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Our  research  group  uses  both  plasma  theory  and  simulation  as  tools  in  order  to  increase  the  under¬ 
standing  of  instabilities,  heating,  transport,  plasma-wall  interaction,  and  large  jmtendals  in  plasma.  We 
perform  plasma  device  computo^  experiments  to  compare  with  analytic  models  and  laboratory  experiments 
in  order  to  accelerate  device  design. 

Our  research  for  1991  has  been  widely  reported,  as  given  by  the  listing  following,  of  Journal  Arti¬ 
cles,  ERL  Reports,  Talks,  and  Poster  Papers. 

Abstracts  are  attached  for  some  of  the  talks. 

Sent  along  with  this  Report  are  reprints  of  Journal  Articles. 

Our  prior  mode  was  to  publish  Quartwly  Progress  Reports;  these  then  became  Semi-Annual  Reports, 
which  end^  in  1988.  In  1989,  we  began  publishing  Annu^  Progress  Reports.  AVhile  QPR’s  were  excel¬ 
lent  exercises  in  reporting,  they  required  an  immense  eRort;  in  today's  research  climate,  such  e^ort  is  not 
available. 

We  trust  that  our  reporting  is  still  useful. 

—  C.  K.  Birdsall 
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PLASMA  THEORY  AND  SIMULATION  GROUP 

Professor  C.K.  Birdsall 

EECS  Department 

University  of  California 

Berkeley,  CA  94720 

U.S.A. 

1.  We  have  developed  the  following  interactive,  real-time,  many  particle  plasma  codes 
for  desktop  computing,  complete  with  graphics: 

•  SPAM:  A  single  particle  mover,  in  specified^,  B  fields;  ruiLS  on  a  PC. 

•  ESI,  XESl:  An  electrostatic  1  dimensional,  periodic  PIC  code.  Can  be  magne¬ 
tized.  Described  in  detail  in  Plasma  Physics  via  Computer  Simulation,  by  Bird¬ 
sall  and  Langdon,  McGraw-Hill  1985  and  Adam  Hilger  1991  (which  has  the  ESI 
disk  with  it).  Input  files  for  projects  in  Chapter  5.  Runs  on  PC’s,  or  on  X-11 
windows  equipped  computers. 

•  EMI:  An  electromagnetic  1  dimensional  code,  periodic  or  bounded.  Not  yet 
fully  debugged. 

2.  We  have  developed  the  following  plasma  device  computer  experiments;  also  for  desk¬ 
top  computing,  complete  with  graphics. 

•  IBC:  Interactive-beam-circuit  code;  a  PIC  traveling  wave  tube,  with  space  charge. 

•  EDPl,  PDCl,  PDSl  (planar,  cylindrical  and  spherical  bounding  electrodes): 
IdSv,  bounded,  magnetized  or  not,  including  electron  and  ion  collisions  with 
neutrals,  external  circuit,  voltage  or  current  sources. 

•  PDP2:  Planar  in  x,  periodic  in  y,  2d3v.  R  -  6  and  R  -  Z  versions  in  early 
development. 

3.  PC  versions  of  SPAM,  ESI  IBC,  PDPl,  PDCl,  PDSl  are  available  at  cost  of  handling 
from: 

Software  Distribution  Office 
Industrial  Liaison  Program 
EECS  Dept.,  Cory  Hall 
University  of  California 
Berkeley,  California  94720 
U.S.A. 

4.  Xgralix(for  X  —  11  windows)  versions  of  SPAM,  EMI,  IBC,  PDPl  and  PDP2  will 
probably  be  available  Spring  1992.  Contact  Prof.  Birdsall  on  these;  do  not  request 
from  ILP. 


September  9,  1991 


PLASMA  DEVICE  COMPUTER  EXPERIMENTS 

or 

PDCX 

Plasma  device  computer  experiments  are  evolving  far  beyond  what  we  have  called  plasma 
simulation  codes  up  until  now^ 

Hence,  let  us  spell  out  what  characterizes  PDCX,  what  is  included  in  PDCX,  cind  some 
of  the  PDCX  methods  and  diagnostics.  This  is  the  purpose  of  this  sheet. 

PDCX  are  characterized  by: 

t 

* 

•  whole  devices:  internal  plasma  and  gas,  external  circuit; 

•  real-time  displays  in  internal  and  external  diagnostics(see  partial  list  at  bottom); 

•  interactive,  in  terms  of  ease  of  viewing  and  re-scaling,  as  well  as  ease  of  setting  initial 
values  of  parameters,  and  in  changing  parameters  any  time  during  an  experimental 
run. 

•  results  in  a  form  to  be  compared  with  laboratory  experiments. 

PDCX  include  the  electromagnetics,  atonaic  physics  and  chemical  reactions  among  charged 
particles(electrons,  ±ions,  particulates)  and  neutral  atoms  and  molecules,  in  three  regions: 

A 

•  surfaces,  walls 

•  sheaths,  edges 

•  bulk  plasma 

PDCX  numerical  methods  involve  a  marriage  among: 

•  PIC,  particle  in  cell,  for  electron,  ion,  neutral  motion,  collective  effects; 

•  MIC,  molecule  in  cell,  for  electron  ion,  neutral  collisions,  short  range  velocity  changes, 
Monte-Carlo; 

•  fluids,  many  kinds; 

•  hybrids  (fluids  plus  particles) 

PDCX  integrations  in  x,v,  t  may  be  explicit  (  to  include  high  frequencies,  short  wave¬ 
lengths)  or  implicit  (  to  keep  only  low  frequencies,  high  wavelengths),  or  be  multi-scaled  in 
I,  V  or  t.. 

PDCX  diagnostics  are  non-invasive,  available  in  i  or  fc  or  t  or  u>,  ajid  are  done  by  species 
(e.g.,  phase  spaces,  n,(i,  t), /,(E,d),(  j*- E),  in  i,t),  plus  the  usual  EM  quantities  and 
collision  statistics,  plus  whatever  the  problem  at  hand  demands. 

Comments  are  most  welcome. 

C.K.  Birdsall,  EECS  Dept.,  University  of  California,  Berkeley,  CA  94720 
September  1991 
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Journal  Articles 


Crystal,  T.L.,  P.C.  Gray,  W.S.  Lawson,  C.K.  Birdsall,  and  S.  Kuhn,  “Trapped  Electron  Effects  on  Time- 
Independent  Negative-Bias  States  of  a  Collisionless  Single-Emitter  Plasma  Device:  Theory  and  Simula¬ 
tion,”  Physics  of  Fluids  B  3(1),  January  1991,  pp.  244. 

Vahedi,  V.,  M.A.  Lieberman,  M.V.  Alves,  J.P.  Verboncoeur,  and  C.K.  Birdsall,  “A  One  Dimensional  Col- 
lisional  Model  for  Plasma  Immersion  Ion  Implantation,”  J.  Appl.  Physics  69(4),  15  February  1991,  pp. 
2008-2014. 

Alves,  M.V.,  M.A.  Lieberman,  V.  Vahedi,  and  C.K.  Birdsall,  “Sheath  Voltage  Ratio  for  Asymmetric  RF 
Discharges,”  /.  Appl.  Physics  69(7)  1  April  1991,  pp.  3823-3829. 

Birdsall,  C.K.,  “Particle-in-Cell  Charged-Particle  Simulations,  Plus  Monte  Carlo  Collisions  with  Neutral 
Atoms,  PIC-MCC,”  IEEE  Trans.  Plasma  Science  19(2),  April  1991,  pp.  65-85  (Invited  paper). 

Procassini,  RJ.,  and  C.K.  Birdsall,  ‘Particle  Simulation  model  of  Transport  in  a  Bounded,  Coulomb  Colli- 
sional  Plasma,”  Phys.  Fluids  B  3(8),  August  1991,  pp.  1876-1891. 

Friedman,  A.,  S.E.  Parker,  S.L.  Ray,  and  C.K.  Birdsall,  “Multi-Scale  Particle-In-Cell  Plasma  Simulation,” 
J.  Comp.  Physics  96,  September  1991,  pp.  54-70. 

Parker,  S.E.,  and  C.K.  Birdsall,  “Numerical  Error  in  Electron  Orbits  with  Large  o)„Sr,”  J.  Comp.  Physics 
97(1),  pp.  91-102,  November  1991. 

Parker,  S.E.,  X.Q.  Xu,  AJ.  Lichtenberg,  and  C.K.  Birdsall,  “Evidence  of  Stochastic  Diffusion  across  a 
Cross-Field  Sheath,”  accepted  for  publication,  Phys.  Rev.  A. 

Parker,  S£.,  A.  Frieoman,  S.L.  Ray,  and  C.K.  Birdsall,  “Bounded  Multi-Scale  Plasma  Simulation:  Appli¬ 
cation  to  Sheath  Problems,”  accepted  by  /.  Comp.  Physics. 

Otani,  NJF.,  J.-S.  Kim,  C.K.  Birdsall,  B.I.  Cohen,  W.  Nevins,  and  N.  Maron,  “Elimination  of  Velocity 
Space  Rings-and-Spokes  Instabilities  in  Magnetized  Electrostatic  Particle  Simulations  of  Plasmas,” 
accepted  by  J.  Coirtp.  Physics  (approx.  June  1992). 

Parker,  S£.,  R.J.  Procassini,  C.K.  Birdsall,  and  B.I.  Cohen,  “A  Suitable  Boundary  Condition  for  Bounded 
Plasma  Simulation  without  Sheath  Resolution,”  accepted  by  7.  Comp.  Physics. 

Verboncoeur,  J.P.,  M.V.  Alves,  and  V.  Vahedi,  “Simultaneous  Potential  and  Circuit  Solution  for  Bounded 
Plasma  Particle  Simulation  Codes,”  accepted  by  J.  Comp.  Physics. 

Berk,  Hi.,  D.D.  Ryutov,  Y.  A.  Tsidulko,  R.H.  Cohen,  and  X.Q.  Xu,  “Electron  Temperature-Gradient  and 
Endloss  Driven  Transport  in  SOL  of  Tokamak  Plasmas,”  submitted  to  Comments  on  Plasma  Physics  and 
Controlled  Nuclear  Fusion. 

Hua,  D.,  X.  Xu,  and  T.K.  Fowler,  “lon-Temperature-Gradient  Modes  in  Non-Circular  Tokamak 
GetMnetry,”  to  be  submitted  by  end  of  year. 

Xu,  X.Q.,  G.  DiPeso,  V.  Vahedi,  and  C.K.  Birdsall,  “Theory  and  Simulation  of  Plasma  Sheath  Waves,”  to 
be  submitted  by  end  of  year. 


Book,  Chapter 


Birdsall,  C.K.,  and  A.B.  Langdon,  Plasma  Physics  via  Computer  Simulation,  Adam-Hilger  edition  (with 
ESI  disk)  1991. 

Birdsall,  C.K.,  and  A.B.  Langdon,  “Particle  Simulation  Techniques,”  in  Computer  Applications  in  Plasma 
Science  and  Engineering,  ed.  A.  T.  Drobot  (Springer-Verlag:  New  York),  pp.  7-41,  1991. 


ERL  Reports 


Parker,  S.E.,  X.Q.  Xu,  AJ.  Lichtenberg,  and  C.K.  Birdsall,  “Evidence  of  Stochastic  Diffusion  across  a 
Cross-Field  Sheath  due  to  Kelvin-Helmholtz  Vortices,”  Memo.  No.  UCB/ERL  M91/79,  September  30, 
1991. 

X.Q.  Xu,  G.  DiPeso,  V.  Vahedi,  and  C.K.  Birdsall,  “Theory  and  Simulation  of  Plasma  Sheath  Waves,” 
Memo.  No.  UCB/ERL  M91/80,  September  30, 1991. 


Conference  Proceedings,  Poster  Papers 


Workshop  on  Edge  Plasma  Physics  for  BPX  andn'ER,Pnn<xuyn,NJ,Janu2ay  15-17, 1991: 

Birdsall,  C.K.,  R.J.  Procassini,  A.  Tarditi,  and  V.  Vahedi,  “1D-3V  Particle  Simulation  of  Tokamak 
Scrape-Off  Layer  and  Divertor  Plasmas” 


International  Sherwood  Fusion  Conference,  Seattie,  WA.  April  22-24, 1991: 

Birdsall,  C.K.,  X.Q.  Xu,  S£.  Parker,  and  AJ.  Lichtenberg,  “Evidence  of  Stochastic  Diffusion  across 
a  Cross-Field  Sheath” 

Hua,  D.,  Fowler,  T.K.,  and  Xu,  X.Q.,  “Gyrokinetic  Particle  Simulation  of  ITG  Modes  in  General 
Toroidal  Geometry” 

Tarditi,  A.,  “2D-Hybrid  Panicle  Model  with  Non-Linear  Electron  Distribution” 

Xu,  X.Q.  G.  DiPeso,  V.  Vahedi,  and  C.K.  Birdsall,  “Theory  and  Simulation  Study  of  Surface 
Waves  in  Bounded  Plasma” 


IEEE  International  Conference  on  Plasma  Science,  Williamsburg,  VA,  June  3-5  1991: 

Tsung,  F.,  J.  Trulsen,  V.  Vahedi,  and  C.K.  Birdsall,  “Simulation  of  Potentials  Created  by  Particu¬ 
lates  in  RF  Discharges:  Residence  at  the  Sheath  Edges” 

Vahedi,  V.,  M.A.  Lieberman,  G.  DiPeso,  C.K.  Birdsall,  T.D.  Rognlien,  J.R.  Hiskes,  and  R.H.  Cohen, 
“An  Atomic  Physics  Model  in  a  Particle-in-Cell  Code  for  Simulation  Plasma  Processing” 

Sixth  International  Corference  on  Emerging  Nuclear  Energy  Systems,  Monterey,  CA,  June  16-21,  1991: 

Avanzini,  P.G.,  and  A.  Tarditi,  “Progress  Towards  a  Neutralized  beam  Experiment  for  a  Colliding- 
beam  Advanced-Fuel  Fusion  Process” 

Tarditi,  A.,  “Multi-Turn  Electron-Ion  Injection  Study  for  Neutralized  High-Density  Beam  Sustaining 
in  a  Closed  Configuration” 


International  Conference  on  Phenomena  in  Ionized  Gases,  Pisa,  Italy,  July  8-12,  1991: 

Tarditi,  A.,  “Particle  Simulation  of  Neutralized  Ion  Bernstein  Waves” 

I4th  Conference  on  Numerical  Simulation  of  Plasmas,  Annapolis,  MD,  September  4-6,  1991: 

Vahedi,  V.,  M.  Surendra,  G.  DiPeso,  and  J.  Verboncoeur,  “Numerical  Methods  for  Simulating  Pro¬ 
cessing  Plasmas” 

Vahedi,  V.,  J.P.  Verboncoeur,  and  C.K.  Birdsall,  “Xgrafix:  An  X-Windows  Environment  for  Real- 
Time  Interactive  Simulations” 


44th  Annual  Gaseous  Electronics  Conference,  Albuquerque,  NM,  October  22-25,  1991 : 

Vahedi,  V.,  P.  Mirrashidi,  B.P.  Wood,  M.A.  Lieberman,  and  C.K.  Birdsall,  “A  Comparison  of  PIC 
Simulation  and  Experimental  Results  in  a  Capacitive  RF  Discharge” 

Lieberman,  M.A.,  V.  Vahedi,  and  R.A.  Stewart,  “An  Analytic  Model  for  the  Ion  Angular  Distribu¬ 
tion  Function  in  a  Highly  Collisional  Sheath” 

Vahedi,  V.,  M.A.  Lieberman,  G.  DiPeso,  and  C.K.  Birdsall,  “A  Bounded  Particle  in  cell  Code  with 
an  Atomic  Physics  Model  for  Simulating  Processing  Plasmas” 

Lieberman,  M.A.,  V.  Vahedi,  and  R.A.  Stewart.  “An  Analytic  Model  of  the  Ion  Angular  Distribu¬ 
tion  Function  in  a  Highly  Collisional  Sheath” 

Mirrashidi.  P.,  B.P.  Wood,  V.  Vahedi,  M.A.  Lieberman,  and  C.K.  Birdsall.  “A  Comparison  of  PIC 
Simulation  and  Experimental  Results  in  a  Capacitive  RF  Discharge” 

Vahedi,  V.,  M.A.  Lieberman,  G.  DiPeso,  C.K.  Birdsall,  T.D.  Rognlien,  J Jl.  Hiskes,  and  R.H.  Cohen, 
“An  Atomic  Physics  Model  in  a  Particle-in-Cell  Code  for  Simulation  Plasma  Processing” 


33rd  Annual  Meeting  of  the  American  Physical  Society,  Division  of  Plasma  Physics,  Tampa,  FL, 
November  4-8, 1991: 

Xu,  X.Q.,  G.  DiPeso,  V,  Vahedi,  and  C.K.  Birdsall,  “Theory  and  Simulation  of  Sheath  Waves  in 
Bounded  Plasmas” 

Vahedi,  V.,  M.A.  Lieberman,  G.  DiPeso,  C.K.  Birdsall.  T.D.  Rognlien,  J  Jl.  Hiskes,  and  R.H.  Cohen, 
“A  Particle  in  Cell  Code  with  an  Atomic  Physics  Model  for  Simulating  Processing  Plasmas” 

Cohen,  R.H.,  and  X.Q.  Xu,  “Scrapeoff-Layer  Instabilities  Driven  by  Temperature  Gradients  and  End 
Loss” 

Tarditi,  A.,  “Hybrid  Particle-Fluid  Simulation  fo  Magnetized  Ion  Plasma-Sheath  Waves” 

Tarditi,  A.,  “Merging-code  approach  for  Realistic  Simulation  of  Plasma  Experiments” 

Theilhaber,  J.,  “Quantum-Molecular- Dynamics  Simulations  of  Liquid  Metals  and  Highly- 
Degenerate  Plasmas”  (Invited  paper) 


1st  Brazilian  Congress  on  Plasma  Physics.  Santos  SP,  Brazil.  December  10-13, 1991: 

Birdsall,  C.K.,  “Particle-in-Cell  Simulations  Plus  Monte-Carlo  Collisions,  PIC-MCC,  for  Partially 
Ionized  Gases,  in  Bounded  Systems”  (Invited  paper) 


Short  Course 


“Plasma  Simulation”  was  taught  to  third  world  professionals  at  the  International  Center  for  Theoretical 
Physics,  Trieste,  Italy,  June  3-13, 1991,  by  C.K.  Birdsall  and  V.  Vahedi. 


Invited  Talks 

C.K.  Birdsall  presented  several  talks  and  live  demonstrations  on  bounded  plasma  computer  experiments: 
Vienna  (Technical  University,  June  18,  1991);  Garching  bei  Miiichen  (Max  Planck  Institute  for  Plasma 
Physics,  July  2, 1991);  Bochum  (University,  July  23, 1991);  Sao  Paolo  (Conference,  December  1991) 


Awards 


1.  C.K.  Birdsall  was  a  Lecturer/Researcher  at  University  of  Innsbruck,  Austria,  February  15  to  July  28, 
1991.  He  presented  a  term-long  course  on  plasma  simulation  and  gave  several  lectures  in  the  local 
plasma  seminar 

2.  At  the  College  of  Engineering  graduation  in  May,  Prof.  Birdsall  was  given  the  Berkeley  Citation,  in 
recognition  of  activity  in  plasma  simulation  and  for  helping  found  the  Energy  and  Resources  Group 
(ERG)  in  1972-1974. 

3.  At  the  14th  International  Conference  on  Numerical  Simulation  of  Plasmas  banquet  on  September  5, 
1991,  in  Annapolis,  MD,  former  students  and  post-doctoral  researchers  of  PTSG  (who  now  number 
almost  SO)  presented  Prof.  Birdsall  a  plaque  (surprise!)  inscribed: 

To  Professor  Charles  K.  (Ned)  Birdsall 

In  warm  appreciation  of  your  many  important  achievements  in  the  sciences  of  electronics, 
plasma  physics,  and  computer  simulation;  of  your  effective  promotion  of  international 
cooperation  in  scieiKe;  and  of  your  contributions  to  the  lives  and  careers  of  the  many  of  us  — 
students,  post-  doctoral  researchers,  and  collaborators  —  who  have  benefited  by  interacting 
with  you,  our  colleague  and  friend. 

CKB  comments* :  This  is  wonderful,  especially  as  it  comes  from  our  simulation  ''family,”  which  has  been 
most  productive  and  very  warm  friends  over  the  past  several  decades.  You  all  have  made  my  33  years  at 
UC  a  very  good  life.  Thank  you! 


*  Yes,  I  »m  now  reiired  from  UC;  however,  while  this  means  no  regular  teaching  schedule.  I  plan  to  continue  in  plasma 
research,  with  PTSG.  for  some  time  to  come.  — CKB 
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Electron  Effects  on  Time-Independent  Negative-Bias  States  of  a  Collision¬ 
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Fluids  B  3(1),  January  1991,  pp.  244 

Vahedi,  V.,  M.A.  Lieberman,  M.V.  Alves,  J.P.  Verboncoeur,  and  C.K.  Birdsall, 
“A  One  Dimensional  Collisional  Model  for  Plasma  Immersion  Ion  Implan¬ 
tation,”  J.  Appl.  Physics  69(4),  15  February  1991,  pp.  2008-2014 

Friedman,  A.,  S.E.  Pailcer,  S.L.  Ray,  and  C.K.  Birdsall,  “Multi-Scale  Particle- 
In-Cell  Plasma  Simulation,”  J.  Comp.  Physics  96,  September  1991,  pp. 
54-70 

Alves,  M.V.,  M.A.  Lieberman,  V.  Vahedi,  and  C.K.  Birdsall,  “Sheath  Voltage 
Ratio  for  Asymmetric  RF  Discharges,”  /.  Appl.  Physics  69(7)  1  April 
1991,  pp.  3823-3829 

Birdsall,  C.K.,  “Particle-in-Cell  Charged-Particle  Simulations,  Plus  Monte  Carlo 
Collisions  with  Neutral  Atoms,  PIC-MCC,”  IEEE  Trans.  Plasma  Science 
19(2),  April  1991,  pp.  65-85  Gnvited  paper) 

Procassini,  R.J..  and  C.K.  Birdsall,  'Particle  Simulation  model  of  Transport  in  a 
Bounded,  Coulomb  Collisional  Plasma,”  Phys.  Fluids  B  3(8),  August 
1991,  pp.  1876-1891 

Parker,  S.E.,  X.Q.  Xu,  A.J.  Lichtenberg,  and  C.K.  Birdsall,  "Evidence  of  Sto¬ 
chastic  Diffusion  across  a  Cross-Field  Sheath  due  to  Kelvin-Helmholtz 
Vortices,”  Memo.  No.  UCB/ERL  M91/79,  September  30,  1991 


III.  ABSTRACTS  OF  1991  TALKS  AND 
POSTERS,  UNPUBLISHED 


Workshop  on  Edge  Plasma  Physics  for  BPX  and  ITER,  Princeton,  NJ,  January 
15-17, 1991  (1  abstract) 

International  Sherwood  Fusion  Cotrference,  Seattle,  WA,  April  22-24,  1991  (4 
abstracts) 

IEEE  International  Conference  on  Plasma  Science,  Williamsburg,  VA,  June  3-5 
1991  (2  abstracts) 

Sixth  International  Conference  on  Emerging  Nuclear  Energy  Systems,  Monterey, 
CA,  June  16-21, 1991  (2  abstracts) 

International  Conference  on  Phenomena  in  Ionized  Gases,  Pisa,  Italy,  July  8-12, 
1991  (1  abstract) 

I4th  Conference  on  Numerical  Simulation  of  Plasmas,  Annapolis,  MD,  Sep¬ 
tember  4-6,  1991  (2  abstracts) 

44th  Annual  Gaseous  Electronics  Coirference,  Albuquerque,  NM,  October  22-25, 
1991  (3  abstracts) 

33rd  Annual  Meeting  of  the  American  Physical  Society,  Division  of  Plasma  Phy¬ 
sics,  Tampa,  FL,  November  4-8, 1991  (6  abstracts) 


WORKSHOP  ON  EDGE  PLASMA  PHYSICS  FOR  BPX  AND  ITER 
l5-n  JANUARY  1991  •  PLASMA  PHYSICS  LABORATORY.  PRINCETON,  NJ 
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3.  SIMULATION  OF  TRANSPORT  IN  A  HIGH-RECYCLING  DIVERTOR  SCRAPE-OFF 
LAYER 
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The  variation  of  the  components  of  (a)  the  electron  and  (b)  the  ion  temperature 
ith  neutral  recycling  intensity  over  the  range  0  <  v'  <  0.25. 


4.  DIVERTOR  PLASMA-SHEATH  REGION  REAL-TIME  SIMULATION 
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5.  aJTUR£  DEVELOPMENTS 
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Theory  and  Simulation  Study  of  Surface  Waves 
in  Bounded  Plasma 

X  Q.  Xu,  G.  DiPeso,  V.  Vahedi,  C.  K.  BirdsaU 
Electronics  Research  Laboratory 
U.  C.  Berkeley,  CA  94720 


Surface  waves  have  been  investigated  analytically  and  with  particle  simulation  for  an 
unmagnetized  2d  plasma  slab  with  periodic  boundary  conditions  in  y  and  bounded  in  x 
with  either  vacuums  (isolated  plasma)  or  conducting  walls.  In  the  vacuum  boundary  case, 
the  particle  simulation  results  are  found  to  agree  reasonably  with  theory  for  high  frequency 
surface  waves.  ^  For  a  plasma  inside  two  absorted  conducting  walls,  simulation  indicate  that 
surface  waves  propagate  along  the  static  sheath-plasma  boundary.  Analytically  treating 
the  sheath  as  a  vacuum  layer,  the  surface  waves  bear  a  resemblance  to  plasma- vacuum 
surface  waves  with  the  vacuum  dielectric  constant  Cq  replaced  by  co  arctan(fcy  A),  where 
A  is  width  of  the  static  sheath.  Nonlinear  interaction  of  bulk  and  surface  waves  in  the 
system  will  be  discussed. 

References 

1.  A.  G.  Sitenko,  Fluctuations  and  Non-linear  Wave  Interactions  in  Plasmas, 
(Pergamon  Press,  Oxford,  1982) 


Evidence  of  Stochastic  Diffusion  across  a  Cross-Field  Sheath 
due  to  Kelvin-Heimholtz  Vortices 

C.K.  Birdsall,  X.Q.  Xu,  S.  E.  Parker^  and  A.J.  Lichtenberg 

Electronics  Research  Laboratory 
University  of  California,  Berkeley,  California  94720 

September  13,  1991 


Abstract 

Our  objective  is  to  identify  the  mechanisms  for  particle  transport  across  a  cross¬ 
field  sheath.  We  present  a  study  of  E  x  B  motion  in  a  vortex  in  which  the  ions  are 
perturbed  by  the  finite  gyroradii  and  electrons  are  perturbed  by  one  or  more  traveling 
waves.  Large  scale  vortices  which  are  the  result  of  a  shear  in  the  E  x  B  drift  velocity 
have  been  observed  in  plasma  simulations  of  the  cross-field  sheath^~^.  Small  scale 
turbulence  is  also  present.  The  vortices  are  the  result  of  the  nonlinear  saturation  of  the 
Kelvin-Heimholtz  instability.  A  vortex  alone  does  not  allow  for  the  observed  electron 
transport  because  the  electron  drift  orbits  simply  circulate.  On  the  other  hand,  the 
ion  motion  can  be  stochastic  from  resonant  interaction  between  the  drift  motion  and 
the  gyromotion,  independent  of  the  background  turbulence.  The  fluctuations  in  the 
ion  density  would  then  give  rise  to  small  amplitude  wave  spectrum.  The  combined 
action  of  the  vortex  fields  and  traveling  wave  fields  on  the  electron  motion  can  agadn 
lead  to  stochasticity.  We  study  these  effects,  showing  that  the  values  of  vortex  fields, 
observed  in  the  simulation,  are  sufficient  to  lead  to  both  ion  and  electron  stochastic 
diffusion.  Furthermore,  the  rate  of  the  the  resulting  diffusion  is  sufficient  to  account  for 
the  diffusion  observed  in  the  simulation. 


Present  address:  Plasma  Physics  Laboratory,  Princeton  University,  Princeton,  New  Jersey  08543. 
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2D-HYBRID  PARTICLE  MODEL  WITH  NON-LINEAR  ELECTRON  DISTRIBUTION 


A.  Tarditi 

Electrical  Engineering  and  Computer  Science,  Electronic  Research  Lab. 
University  of  California  at  Berkeley  -  Berkeley,  CA  94720  (USA) 


ABSTRACT 

A  2D,  hybrid  ( part i c 1 e- i on ,  fluid-electron)  simulation  code 
characterized  by  the  solution  of  the  non-linear  modified  Poisson 
equation,  which  results  assuming  the  Boltzmann  distribution  for  the 
electrons,  is  presented. 

Following  C13,  the  field  solution  is  achieved  through  an  iterative 
procedure.  Anyhow  a  new  scheme  is  considered.  The  potential  is  not 
obtained  by  directly  solving  the  finite  difference  equation  but  via 
the  Green's  function  method. 

The  procedure  begins  with  the  first  guess  for  the  potential.  This  is 
found  through  the  solution  of  the  linearized  modified  Poisson 
equation.  The  Green's  function  for  this  equation,  in  the  2D  case  which 
is  considered,  can  be  found  analytically  in  term  of  .  the  Neumann 
functions  C23,  C31. 

Once  the  potential  corresponding  to  the  linearized  modified  Poisson 
equation  is  known,  the  first  approximation  of  the  electron  (Boltzmann) 
distribution  can  be  calculated.  This  distribution,  plus  t.he  one  given 
by  the  (particle)  ions,  is  considered  as  the  source  term  for  the 
Poisson  equation  (which  now  is  not  "modified"  since  the  fluid  electron 
component  is  taken  into  account  in  the  source  term  itself). 

The  solution  of  this  Poisson  equation  gives  the  second  approximation 
of  the  electric  potential  and  is  still  obtained  via  the  Green’s 
function  method  (as  it  comes  from  the  Coulomb  law,  modified  for  the  2D 
case ,  C  2 1 ) . 

Each  time  step  this  procedure  can  be  iterated  according  to  the  desired 
accuracy.  The  last  iteration  cycle  is  different:  in  fact  the  direct 
solution  for  the  electric  field  can  be  obtained,  without  numerical 
differencing  from  the  potential.  It  is  sufficient  in  this  case  to 
consider  the  electric  field  Green's  functions  (x-  and  y-component  )  for 
the  Poisson  equation  (in  place  of  the  electric  potential  Green's 
f unct ion). 

The  first  results  obtained  with  this  new  code  are  here  presented  and 
compared  with  previous  simulation  runs  based  on  a  linearized  Boltzmann 
distribution  model,  as  in  C23,  C31. 

REFERENCES 

C1]  H.  Okuda,  J.M.  Dawson,  A.T.  Lin  ,  C.C.  Lin,  Phys .  Fluids,  21,  476 

(  1978  ) 

C2I  A.  Tarditi,  Ph . D .  Thesis,  University  of  Genova  (Italy),  1990 
C3]  A.  Tarditi,  submitted  to  the  XX  Int.  Conf.  on  Phenomena  in  Icnized 
Gases,  1991 


Gyrokinetic  Particle  Simulation  of  ITG  Modes  in  General  Toroidal  Geometry 


Daniel  Hua,  T-K.  Fowler 
Department  of  Nuclear  Engineering 
University  of  California 
Berkeley,  CA  94720 

Xueqiao  Xu 

Electronics  Research  Laboratory 
University  of  California 
Berkeley,  CA  94720 


Abstract 


We  have  generalized  the  1-1/2  d  linearized  gyrokinetic  particle 
simulation  code  of  Xu  et.  al.  [1]  to  tokamaks  with  non-circular  corss  section 
(i.e.  elongation  and  triangularity)  to  study  ion  temperature  gradient  modes 
(TTG)  in  the  presence  of  ion-ion  collisions.  We  hope  to  determine  the  role  of 
local  magnetic  shear,  poloidal  beta  and  various  geometric  factors  in  linear 
growth  rate  and  transport  coefficients.  For  comparison  to  Din-D,  we  apply  the 
code  ivith  input  parameters  calcxilated  by  the  MacEquilibrium  Code  of  Haney 
[2]  using  experimental  measurements  of  DIII-D. 

[1]  Xueqiao  Xu,  to  be  published  in  Physics  of  Fluids  £.. 

[2]  Scott  Haney,  Ph.D.  Thesis,  MIT. 


AN  ATOMIC  PHYSICS  MODEL  IN  A  PARTICLE-IN-CELL 
CODE  FOR  SIMULATING  PLASMA  PROCESSING' 


Vahid  Vahedi,  M.  A.  Lieberman 
G.  DiPeso  and  C.  K.  Birdsall 
University  of  California,  Berkeley 
Berkeley,  CA  94720 

T.  D.  Rognlien,  J.  R.  Hiskes,  and  R.  H.  Cohen 
Lawrence  Livermore  National  Laboratory 
Livermore,  CA  94550 


We  are  combining  a  particle-in-cell  (PIC)  model  for  particle  with  a  Monte  Carlo  collision 
(MCC)  scheme  to  model  the  collisions  between  the  charged  and  neutral  particles.  The  MCC 
model  can  also  be  extended  to  model  Coulomb  collisions  between  charged  particles  which  tends 
to  be  significant  at  very  low  temperature  RF  discharges  and  in  ECR  discharges.  We  are  comparing 
the  merits  of  treating  the  neutrals  as  fluids  and/or  as  particles.  These  models  are  incorporated  into 
PDPl^,  a  bounded  one  dimensional  plasma  simulation  code.  As  a  specific  example,  we  consider 
oxygen  RF  discharges,  at  various  neutral  pressures  and  RF  voltages.  The  atomic  physics  model 
for  oxygen  currently  only  includes  the  energy  dependent  processes  of  ionization,  dissociation, 
recombination,  detachment,  and  charge  transfer.  Electrons,  O2*,  O ,  and  O  are  evolved  as  particles. 
We  are  studying  the  modification  of  the  dynamics  of  the  discharge  owing  to  the  presence  of  a 
substantial  concentration  of  negative  ions.  The  atomic  physics  model  can  also  be  used  for  many 
other  types  of  processing  discharges,  such  as  ECR  discharges. 


1.  Work  performed  for  USDOE  by  LLNL  under  contract  W-7405-ENG-48;  a  portion  of  the 
UCB  work  performed  for  NSF  under  grant  ECS-8910827. 

2.  I.  J.  Morey,  V.  Vahedi  and  J.  Verboncoeur,  Particle  Simulation  Code  for  Modeling  Process¬ 
ing  Plasmas,  Bull.  APS.  34:2028,  (1989)  (Abstract);  Codes  available  from  Industrial  Liaison 
Program,  EECS  Dept.,  UC  Berkeley. 
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SIMULATION  OF  POTENTIALS  CREATED 
BY  PARTICULATES  IN  RF  DISCHARGES: 
RESIDENCE  AT  THE  SHEATH  EDGES* 

F.  Tsung,  J.  Trulsen, 

V.  Vahedi,  C.  K.  Birdsall 

University  of  California,  Berkeley 
Berkeley,  CA  94720 


Heavy  particles  may  play  a  role  in  determining  the  average  potentials  experienced  by  ions 
in  RF  discharges,  hence  ion  acceleration  into  targets.  Particulates  or  dust  particles  also  may  play 
a  role  in  many  other  plasmas.  Hence,  it  is  desirable  to  find  where  these  heavy  particles  reside 
(with  respect  to  the  edge  of  the  plasma,  or  sheath)  and  their  effect  of  the  time  average  potential 
which  accelerates  ions  through  the  sheath. 

Using  our  many-particle  PIC-MCC  ld3v  simulation  code  PDPl^,  we  have  been  able  to  show 
that  the  particulates  tend  to  become  charged  negatively,  using  cross  sections  for  electron  and  ion 
attachment  worked  out  here  by  Trulsen,  inspired  by  work  of  R.  Carlile  at  Univ.  Of  Arizona^.  We 
have  placed  one  heavy  particle  at  various  locations  in  the  sheath  and  found  that  special  location 
where  the  time  average  field  at  the  particle  is  zero;  this  is  then  the  residence  of  the  particle,  which 
turns  out  to  be  very  near  the  time  average  sheath  edge.  We  are  now  putting  in  a  large  number  of 
particulates  and  allowing  them  to  move  and  affect  the  potential  across  the  whole  RF  discharge; 
this  is  feasible  only  by  lowering  the  particulate  mass  from  about  10^  argon  ion  mass  to  some 
smaller  values,  still  running  a  long  time  (many  RF  cycles).  We  will  report  on  the  results  at  the 
meeting. 


1.  Work  supported  in  part  by  Univ.  Of  Arizona  Sematech  Ctr.  for  Excellence,  J.  Prince,  Dtr. 

2.  I.  J.  Morey,  V.  Vahedi  and  J.  Verboncoeur,  Particle  Simulation  Code  for  Modeling  Process¬ 
ing  Plasmas.  Bull.  APS.  34:2028,  (1989)  (Abstract);  Codes  available  from  Industrial  Liaison 
Program,  EECS  Dept.,  UC  Berkeley,  CA  94720. 

3.  R.  N.  Nowlin  and  R.  N.  Carlile,  "Electrostatic  Nature  of  Contaminative  Particles  in  a  Semi¬ 
conductor  Processing  Plasma",  ECE  Dept.,  Univ.  Of  Arizona  (Submitted  for  publication) 
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PROGRESS  TOWARDS  A  NEUTRALIZED  BEAM  EXPERIMENT  FOR  A  COLLIDING -BEAM, 

ADVANCED-FUEL  FUSION  PROCESS 


P.G.  Avanaini  and  A. 


Tarditi 


ANSALDO  RICERCHE  -  Corso  Perrone  25  -16161  Genova  (ITALY) 

*Electrical  Engineering  and  Computer  Sciences  -  Electronic  Research  Laboratory 
University  of  California  at  Berkeley  -  Berkeley,  CA  94720  (USA) 


ABSTRACT 

The  problem  of  sustaining  a  colliding 
ion  beam  process  for  advanced-fuel  fusion 
power  generation  is  considered.  In  order  to 
overcome  the  space  charge  limit  on  the  beam 
density,  a  concept  based  on  a  neutralized 
ion-electron  beam  is  introduced. 

A  device  with  fairly  novel  features  is 
described  and  preliminary  design 

considerations  for  a  basic  neutralized  beam 
experiment  are  outlined. 

INTRODUCTION 

An  approach  to  a  fusion  process  based  on 
colliding  beams  has  been  studied  for  several 
years  (1-41. 

The  optimal  reactivity  conditions  for  a  given 
fusion  reaction  are  obtained  when  two 
monoenergetic  counterstreaming  ion  beams 
collide  with  relative  velocity  chosen 
reaching  the  maximum  of  the  fusion  cross 
section. 

The  generation  of  ion  beams  with  energies  of 
hundreds  of  keV,  that  is  the  typical  range 
for  the  reactivity  maxima  of  ver^  energetic 
and  "clean"  reactions  such  as  H-  B  or  D-  He 
15]  is  rather  easy  (while  these  energies  are 
prohibitive  for  confined  plasma  heating  and 
ignition).  A  closed  confinement  geometry  for 
the  colliding  beams  has  to  be  necessarily 
considered  ("open  beam"  configurations  are 
hopeless  for  energy  production  purpouses). 

For  a  reactor  with  reasonable  fusion  power 
density  a  reacting  fuel  density  of  at  least 
10*  m  is  required.  Then,  ion  beams  of  such 
a  density  with  energy  in  the  100  keV  range 
have  to  be  produced  and  confined  while  they 
coll ide. 

At  these  densities,  due  to  the  space  charge 
force  the  beams  must  be  neutralized:  for  this 
purpose  also  a  fast  electron  beam  is 
considered.  For  electrodynamical  stability 
reasons  the  two  ion  beams  run  in  the  same 
direction  at  a  given  relative  velocity  (there 
will  be  a  fast  and  a  slow  beam)  while  the 
electrons  will  be  counterstreaming 
A  “neutralized"  (and  not  “neutral")  beam  will 
then  result  because,  being  charge  balanced. 


it  shows  global  electrical  neutrality  but, 
due  to  the  high  electron-ion  relative 
velocity,  also  low  recombination  (and  elastic 
collision)  rate. 

THE  OPTIMAL  FUSION  PROCESS 

THE  OPTIMAL  REACTIVITY  CONDITIONS 

In  order  to  estimate  the  fusion  power 
production,  the  rate  R«b  (number  of  fusions 
between  the  species  "a"  and  "b"  per  time  and 
volume  unit)  can  be  defined  as 

R  =  n  n  ff  u  (fusions/(m^s) ) 

•b  •  b  r  ab 

where  na,b  are  the  densities  of  the  reacting 
species,  uab  their  relative  velocity  and  <rr 
the  fusion  cross  section.  Here  uab  is  the 
same  for  all  the  particles  (case  of  two 
colliding  monoenergetic  beams). 

The  fusion  power  density  will  then  be 
expressed  as  pf=Rab'Ef  where  Er  is  the  energy 
released  per  fusion. 

The  maximum  rate  is  obtained  when  the 
quantity  oruab  ("reactivity")  is  maximum. 
The  same  definition  holds  when  a 
non-monoenergetic  velocity  distribution  is 
given;  in  this  case  the  reactivity  is  defined 
as  <oTv>  by  averaging  over  the  velocity 
distribution  [4 ) . 

The  curves  of  <rr  versus  particle  relative 
energy  are  well  referred  in  the  literature 
(5).  Then  the  reactivity  vs.  particle  energy 
trend  can  be  readily  found. 

The  reactivity  maxima  for  colliding  beams  and 
thermal  plasmas  are  summarized  in  table  1. 
Here  the  D-T  reaction  has  been  compared  with 
the  most  Important  neutron  free  fusion 
reactions  (characterized  by  high  or  and  Er  at 
the  same  time  [41  ) . 

These  estimates  show  that  the  maximum  ideal 
efficiency  of  the  colliding  beam  approach  is 
greater  than  that  of  the  hot  plasma  roughly 
by  a  factor  of  two. 

In  a  real  device  the  Bremsstrahlung  losses 
must  be  also  taken  into  account:  they  set  a 
limit  to  the  temperature  that  can  be  reached 
in  a  plasma  at  ignition. 


d-CBNtS 


For  a  neutralized  beam  (electron-ion)  beam 
with  axial  velocity  much  greater  than  the 
thermal  velocity  (i.e.  with  small  transversal 
temperature),  the  effect  of  Bremsstrahlung 
radiation  losses  is  less  important. 

COLLIDING  BEAMS  VS.  PLASMAS 


3  3 

Just  pr=10  W/m  (a  very  small  power  density 
for  a  reactor)  the  beam  density  required  for 
each  species  is  n  =(  10^/5. 64- 10  ^)'^^=>4'lo' 
particles/m  .  However  the  stable  confinement 
of  such  a  beam  density  is  unfeasible  due  to 
the  space  charge  limits  and  the  concept  of 
neutralized  beam  has  to  be  introduced. 


The  previous  discussion  about  the 
expected  maximum  fusion  efficiency  must  be 
integrated  by  taking  into  account  the  power 
losses.  The  optimal  conditions  will  be 
reached  when  the  gain  G=pr/pios»  of  the 
process  is  maximum,  where  piosi  is  a  power 
density  that  takes  into  account  the  energy 
absorbed  during  the  process  operation.  Then 
the  real  "optimal"  fusion  cycle  shall  not 
necessarily  correspond  to  the  maximum  fusion 
cross  section  since  the  greater  pr  could  be 
counterbalanced  by  an  increase  in  pio«». 


<0-  V>  1 
r  *  aax 

<r  u| 
r  ' «« X 

E 

f 

(MeV) 

E  <r  u| 
r  r  *M«x 

(Jm^/s) 

<-“b 

3. 5-10‘^^ 
(1  MeV) 

7. 3- 10'^^ 
(700  keV) 

8.64 

10-” 

3-^He 

2.2- 10'^^ 
(300  keV) 

5. 8- 10'^^ 
(420  keV) 

18.  3S 

1.7-  10‘” 

)-T 

8-10‘^^ 
(50  keV) 

210'^‘ 
(100  keV) 

17,59 

5.6'10'^^ 

Table  1.1  -  Reactivity  comparison  for 
colliding  beams  and  thermal  plasmas  in  D-T, 
H-  B  and  D-  He  reactions. 


In  a  thermal  plasma  the  Bremsstrhalung 
radiation  loss  increases  with  electron 
temperature  imposing  then  a  severe  limit  on 
the  energy  gain  of  the  process  and  making 
prohibitive  those  high  temperature  regimes 
required  by  optimal  reactivities  in  advanced 
fuel  reactions  (7J.  In  practice  plasma 
breakeven  will  be  possible  only  at  a  lower 
temperature  than  that  maximizes  the  fusion 
cross  section. 

Then  a  hot  plasma-based  fusion  cycle  will  be 
at  most  a  sub-optimal  process  from  the 
reaction  ideal  efficiency  point  of  view 
(besides  the  problems  concerning  the  heating, 
the  ignition  and  the  confinement  of  plasmas 
at  tens  of  keV  or  more). 

By  trying  to  conceive  a  fusion  process  closer 
to  the  ideal  fusion  cycle  conditions,  a 
colliding  beam  approach  seems  to  be  a 
"natural"  answer. 

As  pointed  out  previously,  the  attainment  of 
ion  beam  energies  of  hundreds  of  keV  does  not 
represent  a  problem. 

The  first  difficulty  arises  from  power 
density  considerations.  For  example,  also  by 
choosing  the  D-T  reaction,  in  order  to  get 


ENERGY  BALANCE  CONSIDERATIONS 

Output  fusion  power  is  produced  as  long 
as  the  Ion  beams  keep  on  running  with 
sufficient  relative  speed. 

The  particles  "burned"  in  the  fusion 
reactions  represent  an  unavoidable  loss  and 
they  shall  be  replaced  in  order  to  keep  the 
density  constant.  Moreover  there  will  be  a 
flux  of  scattered  particles  (ions  and 
electrons)  which,  for  the  high  transversal 
velocity  acquired,  will  tend  to  escape  from 
the  beam.  Some  of  these  particles  will  be 
lost  and  they  have  to  be  replaced,  too. 

Due  to  these  losses  input  beams  supplied  by 
the  injectors  are  also  needed  to  preserve  the 
achieved  regime  conditions  for  a  longer 
time. 

The  system,  in  regime  conditions,  is  then 
continuously  fed  by  two  ion  accelerators  and 
by  an  electron  accelerator  for  replacing  the 
particles  scattered  out  of  the  beam  or 
"burned" . 

The  intrinsic  gain  (i.e.  in  the  ideal, 
lossless  case)  of  the  process  can  be  readily 
expressed  by  considering  the  power  density 
required  to  supply  the  fuel  particles 
"burned"  in  the  fusion  reactions.  If  Ua  and 
Wb  are  the  kinetic  energies  of  the  reacting 
ions,  the  power  density  the  accelerators 
shall  yield  is  pprn=  R»b-(W»+Wb).  where  Rab 
is  the  fusion  rate.  Then  the  intrinsic  gain 
will  be  simply 


p  W  +H 

“^brn  m  b 

By  choosing  Wa-t-UbiUablopt  (the  optimum 
relative  beam  energy  quoted  for  (rru|taax  in 
table  1)  an  intrinsic  gain  of  about  12  for 
the  H-*'b  reaction.  44  for  D-^He  and  176  for 
D-T  is  found. 

This  Is  really  the  maximum  gain  one  could 
expect  since  the  assumption  Wa-^^Wb-Uabjopt 
corresponds  to  the  case  of  counterstreaming 
ion  beams  with  zero  center-of-mass  velocity 
(31 

The  long-range  collisions  produce  a  gradual 
thermalization  and  deteriorate  the 
conf  inement  bui.  they  do  not  throw  the 
particles  out  of  the  beam,  so  particle  losses 
are  duo  directly  only  to  short-range 
col llsions. 

However,  not  all  the  short-range  scattered 
particles  will  be  lost.  This  is  mainly  due  to 
the  fact  that  the  center  of  mass  in  all  these 
collisions  is  moving  (w. r.t.  the  laboratory 
frame)  because  the  ion  beams  are  moving  in 


the  same  direction:  the  scattering 

cross-section  is  defined  for  a  square-angle 
momentum  deviation  in  the  center  of  mass 
frame.  By  means  of  simpie  mechanical 
considerations  [9],  it  is  easy  to  show  that 
when  the  ions  are  moving  fast,  the  scattering 
angle  in  the  laboratory  frame  is  considerably 
reduced. 

BASIC  OPERATION  OF  A  NEUTRALIZED  BEAM  DEVICE 

The  neutralized  beam  process  is  based  on 
two  fundamental  items:  i)  an  injection 

process  that  in  the  start-up  phase  creates  a 
high-density  neutralized  beam  by  accumulating 
particles  from  low-density  electron  and  ion 
Injectors  and  in  the  final  operating  regime 
compensates  the  particles  lost  or  "burned"  in 
the  fusion  reaction;  ii)  a  confinement  system 
which  provides  a  sufficient  beam  stability 
from  the  first  injection  (low-density)  phase 
to  the  final  high-density  pinched  condition. 
The  first  item  has  been  extensively  studied 
and  the  last  results  are  presented  in  (6). 
Direct  supply  of  very  high-density,  100' s  keV 
ion  and  electron  beams  does  not  seem  feasible 
due  technical  and  economic  issues;  it  was 
envisaged  that  a  high-density  neutralized 
beam  can  be  produced  in  a  ring  configuration 
using  low  density,  electrostatic,  ion  and 
electron  accelerators.  The  basic  principle  is 
to  keep  a  sufficient  beam  stability  for  the 
beam  to  allow  a  very  fast  accumulation  of 
particles. 

THE  NEUTRALIZED  BEAM 

THE  ION  AND  ELECTRON  COMPONENTS 

The  basic  idea  is  to  "assemble"  a  high 
density  "neutralized  beam"  by  confining  the 
ion  and  elect, =-ons  together. 

The  electrons  produce  unavoidable  collisions 
and  radiation  losses  leading  towards  the 
thermal izatlon  of  the  whole  system.  However 
the  faster  the  electrons  run  with  respect  to 
the  ions,  the  lower  is  the  electron-ion 
Coulomb  col'lsion  cross  section.  Then  a  very 
fast  electron  neutralization  component  can 
reduce  the  collisional  rate  at  an  acceptable 
level . 

Furthermore  the  high  velocity  electron  beam 
will  have  a  small  energy  spread,  l.e. 
trasversal  temperature,  allowing  a  very  small 
Bremsstrhalung  loss. 

Here  the  ion  energy  is  t/plcally  greater  than 
in  a  discharge-produced  plasma  since  it  is 
provided  independently  to  ion  and  electron 
beams.  The  presence  of  these  high  momentum 
ions  could  have  a  favourable  influence  on 
stability  against  "kink"  or  "sausage"  pinch 
instabi 1 i t ies. 

THE  LOW-DENSITY  COLLISIONLESS  BEHAVIOUR 

The  beam  is  "a  priori"  subjected  to  both 


the  space  charge  fields  produced  by  the  other 
particles  and  to  the  external  confinement 
fields. 

A  simple  estimate,  however,  shows  that  during 
the  injection  phase  the  self-consistent 
fields  'are  negligible  with  respect  to  the 
external  confinement  ones  because  the  low 
beam  densities.  This  means  that  a 
single-particle  treatment  can  be  applied  to 
study  the  neutralized  beam  formation  and 
confinement  in  the  early  phase  of  operation. 
Then,  at  tlje  injection  phase,  the  particle 
trajectories  are  determinate'd  '  nly  by 
external  fields  and  the  electron  and  ion 
beams  constitute  a  non-collisional  system 
which  the  single-particle  treatment  can  be 
applied  to. 

THE  INCREASE  OF  THE  BEAM  DENSITY 

At  the  startup,  low-density  electron  and 
ion  beams  are  injected  into  the  device.  Each 
beam  will  be  constrained  in  a  circular  path 
by  an  external  confinement  and  focusing 
system  [4,  6l. 

In  the  hypothesis  of  perfect  confinement,  if 
L  is  the  length  of  the  beam  orbit,  no  is  the 
beam  density  at  the  accelerator  output  and  vo 
its  (longitudinal)  velocity,  the  density  will 
Increase  linearly  with  time  as 

n(t)  =  n  V  t/L  =  n  t/x 
0  0  0 

being  t=L/vo  the  typical  time  constant. 

For  example,  if  vo=10  ro/s  and  L=1  m,  a 
density  increase  by  three  orders  of  magnitude 
will  take  10‘^  s. 

THE  HIGH-DENSITY  REGIME  CONDITIONS 

The  regime  conditions  for  a  neutralized 
beam  will  be  characterized  by  a  strong 
poloidal  magnetic  field  with  the  related 
pinch  effect. 

If  the  beam  minor  radius  is  much  greater  than 
the  Debye-  lenght  (defined  through  the 
transversal  temperature)  the  neutralized  beam 
is  like  a  neutral,  "soft"  conductor  carrying 
a  strong  current.  In  this  condition  the 
toroidal  curvature  radius  could  be  Imposed  by 
an  external  dipolar  field  directed  as  the 
beam  major  axis. 

The  beam  centrifugal  force  per  volume  unit  is 
m  V  n  +  m  V  n 


where  Rc  is  the  curvature  radius  of  the  beam. 
The  force  due  to  an  external  magnetic  field 
Bo  perpendicular  to  the  beam  plane  is 

F  =  1-lB  , 

m  "^0  , 

where 

1=S- (n  V  q  ♦n  V  q  ) 
e  e  e  1  11 


Is  the  total  current  (S=nRb  is  the  beam 
section)  and  1  is  the  beam  circonference  The 
global  centrifugal  force  on  the  beam  will  be 
Fc=fc-S-1.  Writing  down  the  force  balance 
equation  for  a  one-species  neutralized  beam 
one  gets: 


R  = 

c 


J. 
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m  V  n  m  V  n 

e  e  e  111 


n  V  q  +  n  V  q 
e  e^c  i  i  i 


By  taking  into  account  also  the  "hoop" 
expandlg  force  the  expression  for  Rc  cannot 
be  obtained  anymore  in  a  closed  form  since 


The  solution  "u(t)"  yields  an  estimate  of  the 
effect  of  the  long-range  collisions  on  the 
colliding  beams  velocity  that  is  valid  until 
the  energy  spread  is  low  (since  the 
approximation  of  almost  parallel  velocities). 
A  stability  time  constants  can  be ^defined  as 
the  times  at  which  u(t),  v(t)  or  v  (t)  reach, 
for  example,  the  90%  of  their  initial  value. 
The  short-range  collision  cross-section  is 
simply  (Tir  divided  by  a  factor  41nA. 

A  colllsional  rate  for  short-range  scattering 
can  be  readily  defined  and  taken  into  account 
as  a  loss  term  (the  particles  ^scattered  out 
by  short-range  collisions  could  be  considered 
lost). 


mvn  +  mvn  ul 

e  e  e _ 1  1  1  ^  0  ^ 

nvq+nvq  4jr 

e  e^e  1  II 


if  =  In 


being  "Rb"  the  beam  minor  radius  and  1;  the 
normalized  inductance  per  unit  lenght  and 
li=<H^^>/'H0^(Rb)  (8)  (H0  is  the  poloidal 
field). 

COULOMB  COLLISIOMS 

Coulomb  colllsional  effects  can  be 
studied  by  referring  to  the  long-range  and 
the  short-range  collisions. 

Let  it  be  considered  "test"  charges  q  (mass 
m,  velocity  v,  density  n)  , colliding  wijh 
"field"  charges  q  (mass  m  ,  velocity  v^, 
density  n  ).  The  relative  velocity  is  u=v-v  ; 
the  present  model  holds  until  all  the 
velocities  are  almost  parallel,  then  the 
scalar  notation  will  be  used  for  velocities 
and  momenta.  The  cross-section  for  long-range 
large-angle  collision  is  (8): 

(qq  )^  In  n 

<r  = - 

Ir  .  2  2  4 

4nc  u  u 

where 

In  A  =  b  /b  =  A  47ic  4^u^/(qq  )^ 

aa  X  m I n  D  0 

is  the  Coulomb  logarithm  and  p  the  reduced 
mass.  By  simple  mechanical  considerations  it 
can  be  shown  that  (8) 

^  ''P  ^ 
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THE  "EXTRAP"  HIGH-DENSITY  CONFINEMENT 

When  the  high-density  regime  is  reached 
the  neutralized  beam  looks  like  a  current 
carrying  plasma  in  a  toroidal  z-pinch 
configuration  and  the  confinement  cannot  be 
ensured  anymore  only  by  the  weak-focusing 
effect. 

Lehnert’s  multipolar  plus  z-pinch  "EXTRAP" 
configuration  [10],  seems  to  be  naturally 
suitable  for  the  confinement  of  a 
neutralized,  high-density,  pinched  beam.  Here 
the  EXTRAP  plasma-induced  current  could  be 
replaced  by  the  beam  current. 

It  has  been  experimentally  demonstrated  (10) 
that  EXTRAP  shows  an  excellent  degree  of 
stability  with  respect  to  the  MHD  pinch 
disrupting  modes. 

EXAMPLE 

Some  of  the  most  Important  parameters  for  the 
proposed  process  are  calculated  in  an  example 
case. 

The  D-  He  reaction  was  chosen  and  the  basic 
features  are  listed  below 

Relative  ion  energy  420  keV 

Relative  ion  beam  velocity  8.2-10  m/s 

Energy  released  per  fusion  18.3  MeV 

Fusion  cross  section  7-10”^*  m'^ 

Slow  beam  energy  10  keV  (8-10^  m/s) 

Fast  beam  energy  843  keV  (9-10®  m/s) 

Electron  beam  energy  300  keV  (2.3-10®  ro/s) 
Relativistic  electron  mass  1.58-m 

eO 

1 9  *3 

Assuming  a  regime  ion  density  of  10  m  ,  a 
current  for  the  slow  ion  accelerator  of 
lOOpA,  a  5  mm  minor  radius  beam  and  a  2-m 
long  ring,  the  following  input  parameters  can 
be  found; 


The  time  required  to  reach  the  r^eglme  density 
is  here  2.5  s  and  s  3.8-10  T  external 
magnetic  field  is  needed.  ^  ^ 

The  fusion  power  density  is  1.68-10  W/m  and 
0.78-10*  W/m”^  are  needed  to  replace  the  fuel 
burned.  So  an  intrinsic  gain  of  21.4'  is 
obtained.  The  power  density  required  to 
replace  all  the  short-range  scattered 
particles  is  7-10“  W/m*'’.  then  In  order  to 
get  the  breakeven  at  least  76X  of  these 
particles  must  be  confined  within  the  beam.^ 
The  radiated  power  density  is  about  10  W/m 
(over-estimated) . 

CONCLUSIONS 

A  colliding  beam  fusion  concept  was 
introduced  and  some  fundamental  aspect  of  a 
possible  experimental  device  configuration 
were  discussed. 

The  proposed  approach  is  suitable  for  the 
useof  neutron-free  fusion  reactions.  Besides 
further  theoretical  and  simulation 
investigation,  also  an  experimental  activity 
for  a  basic,  low-cost  feasibility  test  on  the 
neutralized  beam  concept  seems  worthwhile  to 
be  considered. 
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ABSTRACT 

The  problem  of  creating  a  high-density, 
high-current  neutralized  beam  in  a  closed 
configuration  via  a  continuous  particle 
injection  and  accumulation  has  been 
considered  in  relation  to  the  research  on  a 
colliding  beam  fusion  process. 

The  last  results,  oriented  to  the  development 
of  a  basic  neutralized  beam  experiment,  are 
discussed. 

The  theory  of  resistive  wall  injection  of  a 
charged  particle  ring  is  reviewed  and  the 
application  to  the  simultaneous  injection  of 
electron  and  ion  rings  is  introduced. 

INTRODUCTION 

A  colliding-beam  based  fusion  process 
has  been  proposed  and  studied  (1-4);  one  of 
the  major  Issues  for  demonstrating  its 
feasibility  is  to  provide  a  technique  for  an 
efficient  beam  injection  in  a  space-charge 
neutralized  environment:  low-density.  100’ s 
keV-range,  electron  and  ion  beams  have  to  be 
injected  and  trapped  in  a  closed  geometry. 

In  order  to  provide  the  same  equilibrium 
orbit  for  both  the  electron  and  ions,  a 
special  weak-focusing  confinement  system  has 
been  conceived  [31. 

By  maintaining  space-charge  neutrality  and 
through  a  continuous  injection,  a  very  fast 
Increase  in  the  circulating  beam  density 
should  be  achieved,  leading  to  the  formation 
of  a  high-density  neutralized  beam. 

The  injection  process  has  to  be  effective 
from  the  early  stage  of  the  beam  formation 
(low-density)  to  the  final,  high-density 
regime  equilibrium  in  order  to  provide  the 
necessary  compensation  for  the  particles  lost 
(and  "burned"  if  a  fusion  process  were 
implemented) . 

In  order  to  accomplish  these  requirements  a 
"multi-turn",  non-Liouvi Ilian  injection 
technique  based  on  the  axial  motion  damping 
of  electron  and  ion  rings  through  image 
currents  induced  in  a  resistive  cylindrical 
shell  has  been  envisaged.  This  technique  was 
proposed  and  experimentally  tested  in  the 
past  for  ion  and  electron  rings  separately 


(e.g.  [5-8]).  The  present  proposal  is 

intended  to  take  advantage  of  this  previous 
experience  in  a  new  experimental  env,ironment. 

DESCRIPTION  OF  THE  PROCESS 

A  focusing-conf inement  system  for 
holding  electron  and  ion  beams  on  the  same 
stable  orbit  has  been  conceived  [3J.  Then 
low-density  (compared  to  a  fusion  plasma 
density),  counterstreaming,  electron  and  ion 
beams  are  injected  in  a  ring  configuration. 
Each  beam  is  provided  with  a  small  axial 
velocity  in  order  to  miss  the  injector  at  the 
first  turn.  The  ring  motion  generates  image 
currents  in  a  surrounding  resistive 
cylindrical  shell  which  provides  the  required 
damping  (only  in  the  axial  direction,  [6]). 
The  final  condition  is  the  merging  of  the  two 
rings  on  the  same  equilibrium  orbit  resulting 
in  the  formation  of  a  neutralized  beam. 

THE  INJECTION  PROBLEM 

THE  BASIC  THEORY  FOR  THE  INJECTION 

The  fundamentals  of  theory  of  the 
resistive  injection  are  discussed  in  [5]  and 
[6]  for  the  case  of  a  relativistic  electron 
ring.  There  the  solution  of  the  e.m.  problem 
and  the  self-consistent  solution  of  field  and 
motion  equation  for  the  particle  ring  are 
presented.  Those  results  are  here  briefly 
reported. 

The  following  simplifying  assumptions  are 
made  in  order  to  solve  the  problem 
analitically:  i)  the  resistance  of  the 

cylindrical  wall  can  be  approximated  by  the 
surface  resistivity  p,  ii)  the  velocity  along 
’z’  is  considered  non-relativistic  so  that 
the  contribution  of  image  charges  to  the 
retarding  force  is  negligible.  Hi)  the 
current  ring  is  considered  Infinitely  thin. 

In  order  to  design  properly  the  injection 
system  the  equation  of  motion  for  the 
particles  moving  along  "z"  (the  axial 
direction  of  the  cylindrical  shell)  must  be 
solved  with  respect  to  the  beam  (particle 
ring)  parameters. 

In  [6]  the  image  current  density  J(z,t)  on 


the  shell  Is  determined  In  terms  of  the  total 
ring  current  Ir  by  talcing  into  account 
self-consistently  the  motion  of  the  current 
ring  along  "z". 

Due  to  simmetry  reasons  the  only 
non-yaniShlng  components  of  the 
electromagnetic  field  are  Br,  E^,  Bz. 

THE  RETARDING  FORCE  ON  THE  PARTICLE  RING 
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This  last  equation  gives  ( implicitely)  the 
trend  vs.  time  of  the  z-velocity  of  the 
particle  ring  with  respect  the  injection 
parameters. 


The  force  on  each  electron  due  to  the  Br 
component  produced  by  the  image  current  on 
the  resistive' shell  is: 
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In  the  case  the  cylinder  resistance  can  be 
approximated  by  the  surface  resistance  the 
force  depends  only  by  this  last  quantity, 
(i.e.  the  product  of  the  conductivity  <r  and 
the  thickness  “d",  or  o‘-d=d/p  ). 

It  can  be  demonstrated  that,  within  the 
hypothesis  made,  the  image  charge  gives  a 
negligible  contribution  to  the  total  force  on 
the  injected  ring  so  Fh  is  the  only  force 
component  that  will  be  considered  [6]. 
According  to  [5],  for  the  case  of  a  thin 
cylinder  witlj  and  ring  motion  close  to  the 
wall  the  following  expression  for  Fh  can  be 
derived  (written  here  in  MKSA  units): 

c^m  r„  N  2cp/377v  r 
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X  z 

where  me  is  the  electron  mass,  ro  the 
classical  electron  radius,  Ne  the  number  of 
electrons  in  the  ring,  R  the  ring  radius,  “a" 
the  distance  of  this  last  from  the  resistive 
shell  and  377=(p  /c  is  expressed  in  Ohm. 


THE  EQUATION  OF  MOTION  FOR  THE  INJECTED  RING 
The  equation  of  motion  along  "z"  for  one 
electron  (belonging  to  the  particle  ring 
producing  the  image  current  on  the  shell) 
writes  as: 


THE  STOPPING  LENGHT 

An  important  design  quantity  is  the 
"stopping  lenght"  required  to  bring  the  axial 
ring  velocity  to  zero  (reaching  then  the 
equilibrium  position). 

By  dividing  the  eq.  of  motion  by  Vz=dzR/dt  it 
is  obtained: 
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Integrating  by  variable  separation  and 
imposing  vz=0  the  stopping  lenght  Lst=zB-ZRO 
is  found  as 
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(see  eq.  (38)  (51  or  eq.  (21)  in  (6]).  Then: 
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expression  of  Lst  has  to  be  divided  by  v^/c. 
Then,  finally,  by  considering  that  the  ring 
current  is  given  by: 


2  *1/2 

where  yall-(v0/c)  1  in  the  hypothesis 

that  vz«v^. 

Furthermore  if  electric  acceleration  in  the 
z-directlon  is  assumed,  then  the  transverse 
momentum  is  a  constant  of  motion  and  it 
follows: 
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Assuming  yz=l  it  comes 
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By  substituting  the  expression  of  Fh  it  comes: 
2 


N  r  c 

a  0 


377v^/2pc 


y  2  It  R  a  1  ♦  (377v^/2cp) 


This  equation  integrated  by 
separation  yields: 
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while  the  classical  electron  radius  is 
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it  is  obtained 
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It  Is  important  to  notice  that  the  dynamic  of 
the  particle  due  to  the  retarding  force  does 
not  depend  on  its  mass,  *  so  the  stopping 
lenght  for  both  ions  and  electrons  which  are 
in  the  same  current  ring  are  the  same.  The 
quantity  me  in  the  equation  of  Lst  comes. 


throught  expression  of  the  magnetic  field 
produced  by  the  ring  current,  from  having 
Imposed  that  the  electrons  are  confined  with 
a  given  Larmor  radius. 

THE  EXTERNAL  CONFINEMENT  FIELD 

ELECTRONS  AND  IONS  WITH  THE  SANE  LARMOR  RADIUS 

A  vertical  magnetic  field  (with  respect 
the  orbit  plane)  keeps  the  ions  onto  a 
circular  orbit.  A  weak-focuslng  field  profile 
gives  them  the  required  stability. 

When  two  ion  species  are  considered,  if  they 
are  of  different  masses  (like  in  all  the  most 
efficient  advanced-fuel  fusion  reactions  as 
well  as  the  D-f  one)  it  is  possible  to  choose 
the  velocities  of  the  two  species  in  order  to 
both  maximize  the  fusion  cross  section  and 
making  equal  their  Larmor  radii. 

However  the  electrons  will  have  in  general  a 
much  smaller  Larmor  radius  then  the  ions.  The 
only  way  to  keep  them  rotating  on  the  same, 
"reasonable"  size  ion  orbit  (i.e.  in  the 
order  of  the  meter)  is  to  use  relativistic 
electrons. 

By  Imposing  the  equality  of  the  ion  and 
electron  Larmor  radii  it  comes 


where 


smaller  than  the  ion  one  (more  than  one  order 
of  magnitude  with  the  parameters  here 
considered)  and  the  neutralized  beam  could 
not  form  since  the  injection  phase. 

A  radial  electric  field  is  then  Introduced 
(fig.  1).  The  additional  radial  force 
contribution  acts  in  opposite  ways  'for  ions 
and  electrons.  Therefore  by  choosing  proper 
values  for  E  and  B,  the  same  equilibrium 
orbit  for  both  ions  and  electrons  can  be 
imposed. 

In  the  most  general  case,  considering  two 
species  of  charged  particles  with  given 
charge  qe, i,  mass  me,i  and  velocity  ve,i,  the 
force  balance  equations  for  equilibrium 
orbits  with  radius  Re,  i  are: 


then 
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By  Imposing  R=Re=Ri  and  solving  for  E  and  B 
it  is  found  * 
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The  corresponding  electron  kinetic  energy 
will  be  then 


W  =  (7“1)  m  c^ 

e  «0 

For  the  fusion  process  which  it  is  referred 
to  (1-4)  using  relativistic  electrons  with  a 
few  MeV  energy  would  be  unacceptable  from  the 
overall  power  gain  (since  electrons  have  to 
be  supplied  as  well  as  ions  to  compensate  the 
particle  losses). 


THE  "CORRECTED  LARMOR  RADIUS"  CONFIGURATION 

In  order  to  maintain  electrons  and  ions 
on  the  same  orbit  without  being  forced  to  use 
relativistic  electrons  a  special  focusing 
system  was  conceived  [31. 

A  vertical  magnetic  field  B  produces  a 
deflecting  inwards  radial  force,  as  usual. 
The  electron  Larmor  radius  is  however  much 


Fig.  1  -  External  confinement  fields 

For  example,  if  one  assumes  (all  MKSA  units) 
v^  =-10®,  V|  =10®,  R=0.05  for  electrons  and 


protons  respectively,  it  comes  Es2-10^  and 
the  B2l.3-10'^ 

It  must  be  pointed  out  that  the  configuration 
of  E  and  B  here  considered  do  not  produce  any 
ExB  drift  velocity  on  the  particles  because  E 
remains  always  perpendicular  to  the  particle 
trajectory.  Another  way  to  say  this  is 
observing  that  the  particle  guiding  center 
will  remain  always  fixed,  at  the  centre  of 
the  ring. 

WEAK  FOCUSING  WITH  RADIAL  ELECTRIC  FIELD 

The  stability  properties  of  the 
weak-f ocuslng  magnetic  field  profile  remain 
unchanged  with  the  added  electric  component.  A 
charged  particle  (mass  m,  charge  q)  with 
velocity  V  perpendicular  to  a  uniform 
magnetic  field  Bo  will  find  an  equilibrium 
orbit  with  radius  Ro  such  that 

qvB  +mv^/R  =0 
0  0 

The  frame  of  reference  in  fig.  1  is 
considered:  the  equilibrium  o'blt  is  at  r=Ro, 
centered  at  r=0,  the  magnetic  field  has  the 
direction  of  the  2-axls. 

To  acliieve  vertical  focusing,  for  particle 
displacements  in  the  "z"  direction,  the 
magnetic  field  must  decrease  with  r,  that  is 
aB/ar<0,  like  what  is  produced  by  two 
outwards  diverging  poles  (fig  2). 


Fig.  2  -  Magnetic  weak  focusing 

To  ensure  at  the  same  time  the  horizontal 
stability  one  must  impose  that  the 
centrifugal  force  decreases  with  r  faster 
than  the  magnetic  one:  then  a  restoring  force 
component  along  r  will  be  obtained.  This 
condition  is  expressed  as 
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where  (Fa  is  directed  inwards) 


F  =~qvB  ,  F  =mv  /r 


It  comes  then  immediately: 
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By  defining  the  left-hand  side  of  (5)  as  the 
field  index  "n"  and  recalling  the  previously 
stated  condition  aB/ar<0  for  the  vertical 
focusing  (which  is  equivalent  to  n>0)  it 
comes  that  the  stability  of  the  particle 
trajectory  will  be  ensured  if  0<n<l,  that  is 
the  well-known  weak-f ocusing  condition. 

The  same  considerations  will  be  now  applied 
in  the  case  of  the  corrected  Larmor  radius 
configuration. 

The  equilibrium  equation  is 


m  v^/R  =  q  v  B  +  q  E^ 
0  0^0 


If  aE/az=0  the  vertical  focusing  will  not  be 
affected  by  the  presence  of  the  electric 
field  (the  electric  field  is  only  along  the 
radial  direction). 

The  condition  for  the  horizontal  stability 
will  instead  become 
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being  Fe=-qE  the  electric  force  (Fe  must  be 
directed  outwards  for  the  electrons  and 
inwards  for  the  ions). 

Since  it  is  considered  a  radial  electric 
field  produced  by  two  coaxial  cylindrical 
electrodes,  then  E(r)«l/r  which  can  also  be 
written  as  E(r)=EoRo/r,  where  Eo=E(Ro).  So 
3Fe/Sr=qEoRo/r^  and  the  last  equation  becomes 
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Rewriting  this  equation  in' terms  of  the  field 
index  “n"  the  condition  for  the  weak-f ocusing 
will  be  now 
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BASIC  DESIGN  EXAMPLE 

A  calculation  example  Is  reported  for 
giving  an  estimate  of  the  injection 
parameters. 

Electrons  at  300  keV  and  protons  at  100  keV 
are  considered.  Their  velocities^  result 
respectively  v»=2. 3e+08  and  vi=4.  3-10  m/s. 

For  100  keV  protons  electrons  at  13  MeV  are 
required  for  having  the  same  Larmor  radius, 
while  for  protons  at  10  KeV  electrons  at  3 
HeV  are  needed. 

Using  low  energy  protons  could  be  a  way  for 
Implementing  a  neutralized  beam,  low-cost, 
basic  experiment  to  check  the  main  issues  of 
the  proposed  approach  (Injection,  density 
increase,  pinch  effect). 

For  a  electron  beam  injected  at  2- 10®  m/s, 
density  5-10  m~  ,  minor  radius  5  mm  (then 
beam  current  0.12  A)  major  radius  0.2  m, 
distance  from  the  wall  0.01  m,  wall  surface 
resistance  20  n  and  axial  ("z")  velocity  of 
10  m/s,  a  stopping  lenght  of  6  cm  is  found. 
With  a  ring  radius  of  0. 5  m  and  distance  from 
the  wall  5  cm,  a  stopping  lenght  of  0.33  m  Is 
obtained. 

CONCLUSIONS 

A  technique  for  continuous  injection  of 
electron  and  ion  beams  and  their  confinement 
on  the  same  stable  orbit  in  order  to  form  a 
neutralized  beam  was  proposed. 

An  experimental  background  experience  which 
comes  from  previous  experiments  on  electron 
and  ion  beams  separately  can  be  usefully 
applied. 

The  proposed  approach  seems  suitable  for  a 
low-cost  basic  experiment. 
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ABSnWCT 

Tha  Nautralitad  Ion  Barnstain  Uava  INIBU) 
■iaulation  via  an  alactrostatic •  tMo-dinansional 
hybrid  particla-f luid  aodal  is  approachad. 

Siaulations  ara  coaparad  and  eonfiraad  by  alraady 
publishad  axpariaantal  rasults. 

urmBucTiON 

Particla  siaulation  of  aavas  in  plasaasi  dua  to 
tha  intrinsic  aicroacopic  eharactar  of  tha  aodal ■  is 
battar  concarnad  aith  tha  physics  of  tha  partiela-wava 
intaractions  rathar  than  aith  tha  raprasantation  of 
tha  ahola  aava  pattarn  spatial  avolution.  This  is 
particularly  trua  ahan  axplicit  particla  sodals  ara 
usadi  sines  in  tha  spaca  discratization  tha  grid  call 
siza  cannot  ba  such  larqar  than  tha  alactron  Oabya 
langth, 

Tha  particla-fluid  approach  aas  davalopsd  in  ordai  to 
study  plasaa  phanonana  ahich  ara  aainly  affactsd  by 
tha  Ion  dynaaics  Cll.  Tha  high-fraguancy  plasaa 
osrlllatlons  (dua  to  tha  alactrons)  ara  droppad<  tha 
tiaa  stap  siza  is  than  datarainad  froa  tha  ion  aotinn 
tiaa  seals  and  no  anra  froa  tha  alactron  ona. 
Furtharaera  tha  spatial  rasolution  is  noa  linkad  to 
tha  ion  Dabys  langth  in  placa  of  tha  alactron  ona. 
This  aakas  affordabla  problaas  ahich  otharaisa  aould 
raguira  sany  thousands  of  tiaa  staps  and  grid  calls 
pai  linaar  diaansion. 

Mara  ions  ara  siaulatod  as  in  a  convantional  particla 
aodal  nhita  no  coaputar  partielas  ara  usad  for 
alactronst  thay  ara  considarad  lika  a  ‘fluid*  adiich 
aaintains  svaryohara  tha  plasaa  quasi-naiitrality  and 
thair  affact  is  ronsidaiad  only  in  tha  fiSld  solution 
procadura . 

TIC  HYBRID  PMTICXE-FLUIO  CODE  E52MVB 

A  20>  ES  quasi-nautral I  hybrid  (particla  ions, 
fluid  alactrons)  aodal  is  iaplaaantad  in  tha  coda 
F92KfB. 

An  uniforaly  nagnatizad  plasaa  is  siaulatad.  Tha 
aagnatic  fiald  is  along  tha  noraal  diractinn  <zl  aith 
raspact  to  tha  (x-y)  siaulation  plana. 

Tlia  Braan’s  function  'tarhnigua  is  usad  (via  FFT 
convolution)  for  a  diract  solution  of  tha  2-D 
alactrostatic  fiald. 

Tha  fiald  solvar  providas  the  alactrostatic 
salf-consistant  field  produced  by  tha  particles.  Tha 
probleo  is  usually  afforded  by  solving  tha  Poisson 
aguation  (diractlyi  via  Inversion  of  tha  corrasnonding 
sat  of  finite  diffarenca  eguationst  or  through  a 
convolutional  aethod  by  Fourier  transfora  of  tha 
finite  diffaranca  operators,  sea  121)  and  than  by 
nuaarical  differencing  tha  potential  in  order  to 
obtain  tha  electric  field  coaponants. 

Mara  a  Breen’s  function-based  field  solvar  has 
Dnan  davaloped  for  tha  guasi-neutral  hybrid 
particla-fluid  nodal  and  no  diract  rafaranca  to  tha 
Poisson  aguation  is  nada  in  tha  nuaarical  solution 
(tha  potential  nuaarical  gradient  is  also  avoided). 

IMa  electric  fiald  dua  to  a  unitary  2D  charge  in  tha 
origin  Is  considarad  as  tha  Breen’s  function  («  and  y 
roaponants).  Onra  tha  charge  density  distribution  is 
collactad  froa  tha  particla  positions  the  total  fiald 
is  found  after  tha  convolution  batwaan  tha  Breen’s 
function  and  the  charge  density  itself.  Tha 
convolution  is  guickly  perforaen  by  aaans  of  a  2D-FFT. 
According  to  (11  to-  aodel  assuaes  a  Boltzaann 
egnilibriua  for  the  alactrons  with  saall  density 


fluctuations  (a  few  percent)  in  order  to  allOM  the 
linearization  of  tha  Boltzaann  distribution. 

Tha  fiald  solution  proceeds  as  follONs:  first'tha 
aodifiad  Poisson  aguation.  which  takes  into  account 
tha  electron  linearized  Boltzaann  distribution,  is 
salved  analytically  in  two  diaansions  for  a  unitary 
soairce  in  tha  oriqih.  This  gives  tha  Breen’s  function 
for  tha  electrostatic.  Debya-shielded  potential  due  to 
tha  ions.  Than  the  Breen’s  functions  for  the  electric 
field  i-y  coaponants  ara  found  by  analytical 
differencing. 

The  Boltzaann  distribution  for  the  electron  density  n* 
aapanded  at  tha  first  order  writes  as: 

n  >n  axp(eV/k  T  )  Si  n  (l+aU/k  T  ) 

aO  a*  O  Ba 

where  no  is  tha  eguilibriua  electron  density. 

By  considering  a  saall  perturbation  of  the  ion  density 
as  nt>no/Zta6ni  tha  Poisson  agnation  bacoaas 

D^V+k*V— aZ  6n  It 
D  i  k  o 

where  k*-n  a*/k  T  a  ai/\* 

DO  B  a  O  D 

being  Ks  the  electron  Debye  langth. 

Writing  tha  corresponding  hoaogeneous  aguation  for  the 
Breen's  function  in  cylindrical  coordinates  and 
supposing  tha  circular  syaaatry  it  is  found: 


dr 


6  4 
z  V 


dF®: 


♦  k*B 

D  V 


which  is  a  taro-order  Bessel  aiiuation.  It  can  be  shown 
(31  that  tha  Breen’s  function  for  this  aguation.  in 
the  Casa  of  an  isolated  source  .diich  is  considered,  is 


B  (r)  -  V  (k  r)/2n 

V  OP 


whaia  Vo  is  tha  zaro-th  order  Heuaznn  function  (or 
Bessel  functinn  of  tha  second  kind)  in  the  arguaant 
‘kor*. 

In  order  to  find  tha  Breen’s  functions  for  the 
alartric  field  coaponants  it  is  observed  that,  since 
E»  -'TV.  it  coaas 


37  ' 


-  gi  3^  Y  (k  r) 
in  9%  o  D 


Than  it  can  ba  easily  found  that  (31: 


1 

2rr\ 


Y  (?)  - 
1  r 


«y 


1 

2ti\ 

D 


Y,(!)  y/r 


•*>era  Y»  is  tha  first  order  Neuaann  fuiKtion. 

T):a  FFT  convolution  aethod  used  for  tha  fiald 
calculation  assuaas  the  periodicity  of  tha  systao 
•diich  tha  density  and  field  arrays  are  referring  to. 
Tha  coda  is  than  intrinsically  periodic:  the  particles 
passing  over  tha  doaain  border  froa  ona  side  have  to 
ba  Bade  appearing  on  the  other  side  with  tha  saaa 
velocity. 

TW  SInULATION  RESULTS 

It  is  essential  to  define  a  guantitativa 
corraspondenca  batwaan  tha  siaulation  aodel  and  the 
rieal  plaszM  paraaatars.  Tha  coaputar  particles  can  be 
chosen  aith  tha  saae  charga-to-aass  ratio  as  real 
pai  tides,  however,  typically,  tha  values  of  both 
charge  and  aass  are  aany  orders  of  aaqnitude  above  the 


P/S/] 


*  rMl  onr«. 

Th«  definition  of  the  eodcl  por^ootors  is  o  critical 
isouo  for  the  reliability  of  the  sieulation.  A  siople 
alQoritheic  procedure  which  finds  the  input  data  frne 
the  real  plaaaa  specif ications  has  been  used  C31. 

Since  the  external  eagnetic  field  is 
perpendicular  to  the  x~y  siaulation  plane  ion  plasaa 
eaves  propaqate  as  Ion  Bernstein  ec'des.  Horeover  due 
to  the  intrinsic  quasi-neutrality  of  the  eodel 
considered!  only  the  Neutralized  IBM  occurs.  In  the 
case  of  NIBU  in  fact  the  electrons  are  able  to  perfore 
a  sort  of  neutralization  reaching  a  Boltzeann 
equilibriue  with  the  wave  potential  tSl. 

The  propagation  nf  NIBU  nodes  can  be  seen  in  a  narrow 
angular  region  very  close  to  90*  w.r.t.  the  eagnetir 
axis.  A  (i.e.  two-and'half  disensional)  eodel 

would  be  appropriate  for  exploring  the  angular 
dependence  of  the  eode  propagation. 

The  firpt  results  on  the  NIBU  sieulaticn  were  shown  in 
131.  Here  further  calculation  on  the  siaulation  of  a 
real  experiaent  are  reported. 

kixperieental  eeasureaents  on  NIBU  were  recently 
obtained  in  the  (Linear  Hagnetized  Plasaa)  device 
or  HRPP'CPF  in  Lausanne  C61.  The  siaulated  (Argon) 
pljsaa*  as  in  the  experiaent*  has  a  density  of  10 
a  *  Tawl4  eV*  Th«0.1  eV*  with  an  external  aagnetic 
field  of  0.3  T. 


OMEOfl 

ng.  I  -  SKa  eiactric  (veld  eoaciru***  for  moda  • 


r^g.  a  •  VKa  Htsw  dvaparsiOf)  ralaiverc  Okiwulaua''  tm 

eewparad  vitK  IKaory  ond  aaeartwante. 

These  physical  conditions  give  (tipvxuiei^O.  Then*  in 
order  to  obtain  a  sufficient  frequency  resolution  the 
siaulation  has  to  be  extended  over  aany  |/<uc\  but  at 
the  saae  tiae  a  tiae  step  At<l/upi.  (aAt<<l/wev)  has  to 
be  used*  as  the  usual  reqtiireaent  of  explicit  particle 
codes.  Then  a  considerable  nuaber  of  tiae  steps  was 
needed  (20^>. 


The  foraer  siaulation  results  showed  in  C31  are  here 
laproved  by  increasinq  the  nuaber  of  particles  and  the 
spatial  grid  resolution. 

The  code  was  run  on  the  CRAY-2  with  128  •'lEG  grid 
points  and  6SS36  (*256  >  particles. 

In  fig.  1  an  exaaple  for  the  electric  field  spectrum 
(for  the  3*'^  aode)  is  shown.  Up  to  eight  aodes  were 
recorded  with  k*a/c!pc^» 

The  dispersion  relation  is  obtained  as  in  fig.  2.  The 
experiaental  results  and  theoretical  curves  (both  froa 
C61)  are  reported  on  the  saae  plot  for  coaparison.  The 
frequency  resolution  was  A(.>3suMi./4.St  2066  tiae  steps 
with  At*2  10  having  been  used.  A  very  good  agreeaent 
IS  found  for  the  first  two  branches.  The  dif/erences 
in  the  third  one  are  attributed  to  the  still  not 
suffiriently  large  scale  of  the  siaulation  considered 
so  far. 

A  aaxwellian  plasaa  was  initialized  through  a  "quiet 
start"  technique  and  a  good  stability  of  the  theraal 
equiltbriua  was  obtained  as  showed  by  the  final 
distribution  function  in  fig.  3  (after  about  600 
plasaa  peiiods).  « 


V  MCOuiJS 


r^g.  9  •  tor  diitrvt>wti.or  ot  tK«,  ord  of  run 


OMCLUStONS 

A  2D,  ES  quasi-neutral  , hybrid  porticle-f luid 
code.  Mas  applied  to  the  siwilation  of  NiBw. 

The  hybrid  code  exploits  a  new,  rather  sophisticated, 
field  solver  for  a  linearized,  quasi-neutral,  plasaa 
aortel.  The  extension  to  the  non-linear  case  could  be 
Made  straightforwardly  by  successive  field  solving 
withiu  each  tiae  step.  Uc*rk  m  thix. 
direction  is  in  progress. 

Techniques  for  defining  correlations  between  the 
particle  eodel  and  the  real  pla'^ma  were  exploited. 

The  character  of  the  siaulation  results  is  still  that 
of  an  approach  phase.  Anyhow  it  is  believed  that  a 
contribute  to  the  ieproveaent  of  the  particle 
siaulation  aethod  reliability  and  to  the  progress 
towards  a  closer  link  with  plasM  experiaents  has  been 
given. 
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RF  glow  discharges  and  other  processing  plasmas  are  used  extensively  in  the  microelec¬ 
tronics  industry.  Self-consistent  fluid  equations  have  been  used  recently  to  study  the  structural 
features  of  RF  and  DC  glows*'^.  However,  since  these  discharges  are  inherently  complex  and  the 
particle  distributions  are  non-Maxwellian,  there  has  been  a  considerable  effort  in  developing 
self-consistent  kinetic  models  without  making  any  assumptions  on  the  distribution  functions'^. 


In  order  to  use  particle-in-cell  simulation  codes  for  modeling  collisional  plasmas  and  self- 
sustained  discharges  it  is  necessary  to  include  interactions  between  charged  and  neutral  particles. 
Monte  Carlo  methods  have  been  used  extensively  in  swarm  simulations*'*^.  In  many  Monte  Carlo 
schemes,  the  time  (or  distance)  between  collisions  for  each  particle  is  calculated  from  a  random 
number.  This  allows  for  efficient  algorithms,  especially  when  the  null  collision  method  is  used*^. 
This  technique  is  however,  not  compatible  with  PIC  simulations,  since  all  particle  trajectories  are 
integrated  simultaneously  in  time.  Thus,  the  collision  probability  for  the  ith  charged  particle  is 
calculated,  based  on  the  distance  As,  =  v,Ar  travelal  in  each  time  step  At,  to  be 

F,  =  1  -exp(-As,ar(£,)n) 

where  is  the  total  collision  cross  section,  £,  is  the  kinetic  energy  of  the  particles  and  n  is  the 

neutral  density.  A  collision  takes  place  if  a  uniformly  distributed  random  number  on  the  interval 
[0, 1]  is  less  than  F,.  The  null  collision  method  can  be  incorporated  into  the  collision  model  by 
picking  a  constant  collision  frequency  v'  such  that,** 

v'^/IVOr 

Thus,  the  computational  cost  of  calculating  F,  can  be  avoided.  The  fraction  of  the  total 
number  of  particles  (chosen  at  random)  in  the  simulation  that  experience  collisions  is  given  by 

F  =  1  -exp(-v'A/) 


The  collision  is  assumed  to  take  place  at  the  current  position  of  the  particle.  It  should  be 
noted  that  the  choice  of  At  will  affect  the  accuracy  of  the  collision  model.  For  instance.  As,  should 
not  be  much  larger  than  simulation  length  scales  of  interest  (e.g.  grid  spacing,  X-o,)  and  As,ar(£,  )/i 
should  be  about  0.1  or  less*’.  Once  a  collision  occurs,  the  type  of  collision,  the  energy  of  the 
ejected  electron  (for  an  ionizing  collision)  and  the  direction  (s)  are  determined  with  new  random 


numbers.  These  quantities  are  related  back  to  the  system  coordinate  axes.  The  procedure  for 
electron-neutral  collisions  is  describe  in  detail  by  Boeuf  and  Marode'”,  and  by  Thompson  et  al^* 
for  ion-neutral  collisions.  Expressions  for  differential  cross  sections  that  are  analytically  integrable 
are  useful  as  the  computational  cost  of  determining  scattering  angles  and  energy  redistribution  in 
ionizing  collisions  is  minimized^’ 

A  Monte  Carlo  collision  handler  as  described  above,  including  the  null  collision  method, 
has  been  developed  as  an  addition  to  the  PIC  scheme  as  shown  in  Fig.  1 .  The  full  three  dimensional 
character  of  a  collision  is  modeled  with  three  velocity  components.  The  neutrals  are  assumed  to 
be  uniformly  distributed  between  the  boundaries  with  a  constant  density  and  a  Maxwellian  profile. 
The  model  is  still  valid  if  the  neutral  density  is  a  weak  function  of  position  and  time  (small  variations 
across  the  mean  free  path  and  collision  times).  This  scheme  can  also  be  extended  to  model  Coulomb 
collisions  between  charged  particles. 


FIG.  1 .  The  flow  chart  for  an  explicit  PIC  scheme  with  the  addition  of  the  collision  handler,  called 
PIC-MCC*. 

RF  discharge  modeling  displays  many  physical  time  scales,  e.g.  <  0)^,.  With  a  PIC 

model  including  an  electrostatic  response,  the  highest  frequency  that  must  be  resolved  by  the 
explicit  numerical  methods  used  to  solve  the  particle  and  field  equations  is  co^,.  If  (Op^A/  >  1, 
numerical  instabilities  can  occur  for  explicit  methods**.  To  observe  the  physics  of  interest,  one 
needs  to  resolve  th  RF  timescale  only,  and  therefore  much  computing  time  is  wasted  resolving 
the  plasma  oscillation  time  scale. 

Implicit  particle  simulation*’  has  been  developed  to  relax  the  numerical  stability  constraint. 
We  will  now  briefly  review  implicit  particle  simulation  techniques.  Implicit  particle  movers 
advance  the  postion  of  the  <th  super  particle  by  the  equation 


where  x  is  the  portion  of  the  position  advanced,  dependent  on  quantities  known  at  present  and 
previous  time  levels,  a'  =  and  P  depends  on  the  particular  implicit  scheme.  The  field 

at  a  particle  location  is  interpolated  from  the  field  known  at  the  grid  in  space.  Note  that  the  particle 
location  at  the  future  /i+l  time  step  depends  on  the  field  at  that  time  step,  but  the  field  at  the  future 
time  step  depends  on  the  particle  location  at  that  time  step  through  the  Poisson  equation. 

One  way  to  get  around  this  problem  is  to  linearize  the  locations  x,"  ‘  about  the  locations 
jc"  *  ‘  in  the  superparticle-to-grid  weighting  equations.  Then 


where  p"^ '  is  determined  by  weighting  superparticles  at  jc,"*'  to  the  grid.  The  minus  sign  is  from 
the  functional  dependence  of  the  particle  positions  in  the  superparticle-to-grid  weighting  equations. 
Strictly  speaking,  5x"  *  ‘  =  Jt"  *  *  -  jt "  * '  is  an  individual  quantity  for  each  superparticle  and  is  given 
by  the  advancing  equation.  Instead,  =  5x"*‘  =  i.e.,  the  perturbations  are  taken  to 

be  grid  quantities  while  maintaining  the  form  as  given  by  the  advancing  equation.  Inserting  the 
above  form  of  the  perturbation  into  the  equations  for  the  density  and  combining  with  the  Poisson 
equation  gives  the  numerically  implicit  Poisson  equation 


where  a  =  a'/c,,. 

The  equation  is  solved  on  a  spatial  grid.  Simple  boundary  conditions  for  the  RF  discharge 
are  a  zero  potential  at  the  left  wall  and  an  RF  source  voltage  at  the  right  wall.  The  electric  field 
is  found  at  the  interior  points  by  central  differencing  the  potential.  The  electric  field  at  the  walls 
is  given  by  a  numerically  implicit  Gauss’  law  which  is  the  integral  form  of  the  numerically  implicit 
Poisson  equation.  The  electric  field,  at  the  left  wall  for  example,  is  then  derived  by  taking  a 
Gaussian  pill  box  about  the  wall, 

[(1  -  ap)E]^ ,  =  ao/e„  poAx/2£o 


The  enclosed  charge  (RHS)  includes  the  wall  surface  charge  density  Oq  =  where  the  0  subscript 

indicates  the  left  wall  quantity  and  j= 1/2  indicates  an  evaluation  between  the  0  and  1  grid  points. 
It  is  then  possible  to  solve  for  the  electric  field  on  the  left  wall.  A  similar  proceedure  is  used  for 
the  right  wall.  The  equations  may  be  generalized  to  included  more  complicated  boundary  con¬ 
ditions  including  external  circuit  elements. 

Many  accuracy  constraints  still  remain.  One  important  constraint  is  that  the  fastest  particle 
species  should  resolve  spatial  gradients  in  the  field,  i.e.,  <  1  where  s  is  the  gradient  length. 

Another  important  constraint  is  that  all  particles  should  sample  the  field  on  the  grid  in  a  continuous 
manner  over  a  single  time  step.  This  gives  v^^AiIAk  <  1.  A  problem  with  implicit  methods  is 
excessive  numerical  cooling  which  is  due  to  poor  sampling  of  fast  particles  in  simulations  with 
large  time  steps.  Resolving  fast  particles  is  particularly  important  in  RF  discharges  because  it  is 


the  fast  electrons  which  maintain  the  discharge  through  ionization  collisions  with  the  neutrals.  A 
possible  way  out  of  this  problem  is  to  do  multi-scale“  simulations.  That  is,  the  few  fast  electrons 
that  maintain  the  discharge  may  be  pushed  with  a  small  time  step  while  the  remaining  slow  particles, 
essentially  residing  in  the  bulk  plasma,  may  be  pushed  with  a  large  implicit  time  step. 
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We  have  developed  a  real-time  user  interface  environment,  XGRAFIX,  for  simulations 
running  under  X-Windows.  XGRAFIX  is  written  in  C  in  an  object  oriented  style,  and  since  it 
uses  only  the  lower  level  X-Windows  function-calls,  it  can  be  compiled  with  any  superset  of 
X-Windows,  e.g.  Motif,  and  is  compatible  with  many  systems.  XGRAFIX  allows  the  .user  to 
display  multiple  diagnostics  and  view  them  as  they  evolve  in  time.  Like  most  other  X  environments, 
XGRAFIX  provides  keyboard  and  mouse  suppoits.  The  simulation  codes  are  structured  as  shown 
in  Fig.  1. 


Figure  1 .  Schematic  representation  of  the  interaction  between  XGRAFIX  and  the  physics 

kernel. 


The  physics  kernel  is  portable  to  any  machine  supporting  standard  C.  The  INTT  module 
scans  the  input  file  containing  the  physical  parameters  of  the  problem  and  initializes  the  diagnostic 
windows.  INTT  also  sets  up  memory  for  array  storage.  The  environment  provides  hooks  for  the 
physics  kernel  to  run  continuously  (there  is  no  time  limit  -  the  code  can  run  indefinitely)  or  step 
through  individual  timesteps.  The  windows  are  refreshed  each  timestep,  and  all  user  requests  are 
processed  by  the  XGRAFIX  Event  Manager.  When  the  simulation  is  in  the  running  state,  the 


Event  Manager  is  invoked  at  each  time  step  to  process  the  events;  if  not  running,  the  Event  Manager 
is  constantly  invoked.  The  events  include  moving,  resizing,  and  iconifying  windows  as  well  as 
mouse  button  clicks  and  keyboard  inputs.  XGRAFIX  also  supports  PostScript  output  of  the  plots. 

The  environment  currently  offers  three  types  of  plots;  linear,  semi-log,  and  scattered. 
XGRAFIX  provides  the  user  with  menus,  dialogboxes,  and  smart  windows.  Each  window  has 
(currently)  four  standard  buttons  for  rescaling  the  graph,  viewing  traces  of  plots,  PostScript  output, 
and  a  crosshair  for  measurements.  The  menus,  dialogboxes,  and  the  mouse  make  the  environment 
especially  user-ftiendly. 

XGRAFIX  is  being  currently  rewritten  in  C++  which  offers  pre-defined  classes  and  hier¬ 
archies  used  in  the  object-oriented  style  for  such  objects  as  buttons,  menus,  windows,  etc.  This 
modification  should  make  it  easier  to  handle  future  additions  to  XGRAFIX.  We  are  also  adding 
optimized  three  dimensional  plotting  routines  to  the  environment  which  will  make  XGRAFIX 
even  more  useful  for  2D  and  3D  simulations. 

We  are  presently  running  three  of  our  codes,  ESI  (Electrostatic  1  Dimensional  periodic 
plasma  simulation)^,  PDPl  (Plasma  Device,  Planar  1  Dimensional)*  and  PDP2  (Plasma  Device, 
Planar  2  Dimensional)  in  XGRAFIX. 

This  work  will  be  presented  at  the  74**  International  Conference  on  the  Numerical  Simulation 
of  Plasmas. 
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mean  ficee  paths  for  ion-neutral  scattering  and  charge-exchange 
collisions,  respectively.  With  this  assumption,  we  consider  the 
angular  distribution  to  arise  mainly  from  ions  that  strike  the  electrode 
after  undergoing  only  one  scattering  collision  following  the  last 
charge-exchange  collision. 
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collisions  between  the  charged  and  neutral  particles.  The  MCC  model 
can  also  be  extended  to  model  Coulomb  collisions  between  charged 
particles  which  tends  to  be  significant  at  very  low  temperature 
discharges.  These  models  are  incorporated  into  PDPl,  a  bounded 
one  dimensional  plasma  simulation  code.  As  a  specific  example,  we 
consider  oxygen  RF  discharges  at  various  neutral  pressiues  and  RF 
voltages.  Electrons,  O/,  O',  and  O  are  evolved  as  partilces.  These 
models  can  be  used  to  model  other  processing  discharges. 
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Bi&osall,  UniversUy  of  CaHfomia,  Berhdey — ^Sheath  waves  have  been 
investigated  analytically  and  with  partick  itiTniila.tirwi  {qx  an  nnmagne- 
tized  two  dimensional  plasma  slab  with  periodic  bonndaiy  r-rniHitjmia 
in  y  and  conducting  wadis  in  x  Analytically  treating  the  sheath  as  a 
vaicnnm  layer,  the  sheath  wave  bears  a  resemblance  to  plaOTna  vacnnm 
snziace  waves.  The  simnlations  are  in  good  agreement  with  the  theory 
for  both  bnlk  Bohm  Gross  waves  and  ^ge  sheath  waves.  We  have  also 
simulated  a  magnetized  plasma  in  both  the  pure  (PIBW)  and  neutral¬ 
ized  (NIBW)  ion  Bernstein  wave  ze^mes  to  look  for  sheath  waves  in 
these  cases.  For  PIBW  and  NIBW,  the  ions  are  fully  magnetized  while 
for  PIBW,  the  electrons  Sixe  treat^  as  a  background  and  for  NIBW, 
they  are  treated  in  the  drift  kinetic  approximation.  Ultimately,  we  want 
to  do  magnetized  simulations  to  gain  understanding  of  the  impurity  and 
edge  heating  problem  for  ICRF  experiments. 
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Lawrence  Livermore  NaHonal  Laboratory — We  are  combining  a  particle 
in  cell  (PIC)  model  for  particle  and  fidd  dynamics  with  a  Monte  Carlo 
(MCC)  scheme  to  model  the  atomic  physics  of  particle  collisions  with  a 
background  neutral  gas.  These  modds  are  incorporated  into  PDPl,  a 
one  dimensional  bounded  plasma  simulation  code.  As  a  specific  example, 
we  consider  oxygen  RF  discharges  at  verious  neutral  gas  pressures  and 
RF  voltages.  Electrons,  O^,  O",  and  O  are  evdved  as  partides.  These 
models  can  be  used  for  other  processing  discharges.  Due  to  the  varying 
time  scales  of  dectric  fidd  and  coUisional  dynamics,  discharge  equilibri¬ 
um  is  difficult  to  reach  in  the  simulations.  We  will  discuss  methods  to 
acheive  equilibrium  more  rapidly. 
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Scrapeoff-Layer  Instabilities  Driven  by  Tempera¬ 
ture  Gradients  and  End  Loss  T  R.  H.  Cohen,  LLNL;  X.  Q.  Xu, 
U.C.  Berkeley — We  have  performed  a  kinetic  analysis  of  an  instability^ 
driven  by  electron  temperature  gradients  in  the  presence  of  end  loss,  and 
examined  its  applicability  to  tokamak  scrapeoff  layers.  This  instability 
is  strong  enough  that  it  could  set  the  width  of  the  scrapeoff  layer.  In 
addition  to  kinetic  effects,  our  analysis  adds  secondary  electrons,  recy¬ 
cling,  energy  endloss  and  (where  appropriate)  electron  Landau  damping, 
and  additional  finite-gyroradius  terms.  The  dispersion  relation  is  un¬ 
changed  significantly  by  the  transition  &om  a  collision-dominated  fluid 
regime  to  a  more  collisionless  kinetic  regime,  so  long  as  the  ordering  ^ 
Vti/L||  <<  u;  <<  vte/X||  and  a  lower  bound  on  collisionality  are  satisfied, 
but  the  energy  endloss  introduces  a  threshold  in  L\\/Lt-  The  atomic- 
physics  corrections  reduce  growth  rates  and  raises  the  threshold.  We 
re-analyze  the  mode  with  alternate  orderings.  For  realistic  parameters 
the  ordering  expansions  are  marginally  justified;  furthermore,  axial  gra¬ 
dients  in  equilibrium  quantities  are  significant.  For  these  reasons,  and 
to  begin  to  assess  nonlinear  effects,  we  have  developed  a  2D  gyrokinetic 
simulation  model.  We  report  analytic  and  simulation  results  for  Dili 
and  ITER  parameters. 
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Hybrid  partide-fluid  simulation  of  magnetized  ion  plasma- 
sheath  waves.  A.  Tarditi,  EECS  University  of  California  at  Berkeley, 
Berkeley,  CA  (USA).  A  2D  electrostatic  particle-ion,  fluid-electron  code 
has  been  developed  for  studying  ion  waves  in  a  magnetized  ion-sheath 
region  close  to  an  absorbing  conducting  wall.  A  semi-infinite  plasma 
slab  model  is  considered,  x-bounded  and  y-peiiodic.  In  a  hybrid  model 
the  particles  (ions)  are  advanced  on  the  ion  time-scale  while  a  fluid-l^e  , 
electron  component  is  simulated  by  providing  an  analytic  expression 
for  a  non-uniform  density  which  t^es  into  account  the  overall  effect 
of  the  guiding  center  electron  diffusion.  Then  in  the  sheath  region  a 
non-Boltzmann  electron  density  profile  (function  of  the  plasma 
potential),  which  is  eventually  merged  with  the  Boltzmann  distribution 
in  the  "bulk"  plasma  region,  is  considered.  An  iterative  field  solver  for 
a  non-linear  Poisson  equation  has  been  implemented  in  the  code,  so 
virtually  any  expression  of  electron  density  profile  vs.  potential  can  be 
considered. 

Work  supported  by  NATO  Advanced  Fellowship  Program,  US  DOE  contract  DE- 
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Merging-code  approach  for  realistic  simulation  of  plasma 
experiments.  A.  Tarditi,  EECS  University  of  California  at  Berkeley, 
Berkeley,  CA  (USA).  A  merging-code  approach  for  complex,  space- 
time  multiscale  simulation  of  plasma  experimental  devices  is  proposed. 
This  study  is  particularly  oriented  to  the  Numerical  Tokamak 
experiment  (NTX)  project  Different  codes,  operating  on  different 
space  and/or  time  scales  as  well  as  dealing  with  different  sets  of  physical  * 

parameters,  run  at  the  same  time  while  sharing  and  exchanging 
interactively  their  I/O  data  fluxes.  The  information  interchange  is 
p^ormed  through  physical  quantities  (i.e.  density  profiles,  current 
distribution,  etc.)  ratiier  than  through  numerical  quantities  (late  com¬ 
puter  particles)  which  are  code-dependent 
Physics  and  Computer  Science  issues  relevant  to  the  development  of  a 
mititi-task  integrated  simulation  are  discussed.  The  features  of  a 
versatile  software  environment  (shell),  based  on  a  merged-code  stan¬ 
dard  shell  communication  protocol  and  provided  with  a  modular 
structure  which  allows  the  use  of  even  previously  written  codes,  are 
described.  A  simple  simulation  example  is  also  referenced. 
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lon-Temperature-Gradient  Modes  in  Non-Circulax 
:  Geometry.*  D.  Hua,  X.  Xu,  and'T.K.  Fowler,  UC  Berkeley. 

■  —  A  l|-D  linearized  gyrokinetic  code  employing  a  5f  particle 

;  algorithm  for  ion  temperature  gradient  mode  cadculations  is  ex-  . 
tended  to  non-circular  cross-section  to  study  scaling  with  the 
elongation  k.  The  growth  rate,  wave  numbers  amd  stability 
thresholds  are  calculated  for  re  ranging  from  1  to  2,  for  a  fiat 
density  profile  characteristic  of  H-mode  operation.  Calibration  : 
with  a  hot-ion  H-mode  shot  in  DIII-D  gives  fair  agreement  be-  . 

;  tween  mixing  length  estimates  of  Xt  derived  from  the  theory  and 
;  experimentally-derived  values.  Agreement  requires  taking  into 

■  account  Tg  Ti  in  the  core,  which  reduces  Xi  neajr  the  axis. 
However,  the  theory  fails  to  account  for  the  strong  dependence 

'  on  rc  at  fixed  q  implicit  in  empirical  scaling  laws  for  the  energy 
confinement  time. 
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Abstract 

In  recent  years  particle-in-cell  techniques  with  Monte  Carlo  collisions  have  shown  to  be  of 
great  use  in  modeling  discharges  and  processing  plasmas  [1]  [2]  [3]  [4],  A  better  understanding 
of  these  plasmas,  e.g.  RF  driven  discharges,  has  enabled  us  to  device  a  scheme  to  optimize  the 
modeling  code.  It  will  be  shown  how  direct  implidt  particle  simulation  [5]  and  multi-time  scale 
scheme  [6]  are  incorporated  into  the  bounded  one  dimensional  plasma  simulation  code  PDPl  in 
order  to  relax  the  UptAt  time  constraint. 


1  Introduction 

In  planar  RF  discharges,  a  plasma  is  bounded  between  two  parallel  plates  and  is  driven  by  an 
external  source  as  shown  in  Fig.  1.  The  electrons  responding  to  the  instantaneous  applied  held  gain 
energy  by  colliding  with  the  moving  sheath,  and  the  electron  distribution  develops  a  high  energy 
tml.  A  typical  electron  energy  distribution  function  is  shown  in  Fig.  2.  Although  the  population 
of  the  high  energy  tml  is  down  by  several  order  of  magnitude  from  the  bulk  plasma  population,  it 
is  the  electrons  on  the  tail  of  the  distribution  which  overcome  the  ionization  threshold  to  keep  the 
discharge  alive  through  ionization  collisions  with  the  neutral  particles.  In  modeling  discharges,  one 
must  pay  s  careful  attention  to  resolve  the  orbit  of  the  high  energy  electrons  accurately. 

With  a  particle  simulation  model  based  on  an  electrostatic  plasma  response,  i.e.  particle  equa¬ 
tions  of  motion  coupled  to  the  Poisson  equation,  the  highest  frequency  that  must  be  resolved  by 
the  explicit  numerical  methods  used  to  solve  the  particle  and  field  equations  is  Up^.  That  is,  if 
UptAi  >  1,  numerical  instabilities  can  occur  [7].  Direct  implicit  particle  simulation  [5]  relaxes  the 
UptAi  time  constrrunt.  However,  The  temptation  to  use  a  larger  At  is  foiled  because  the  accuracy 
condition  vAt/Ax  <  0(1)  would  be  violated  for  the  fast  electrons.  The  poor  sampling  of  the  fast 
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electron  orbits  may  lead  to  numerical  cooling  of  the  high  energy  tail  of  the  distribution  making  it 
impossible  for  the  electrons  in  the  model  to  reach  ionization  threshold. 

A  multi-time  scale  scheme  is  deviced  to  allow  the  fast  electrons  to  move  with  a  small  time  step 
while  the  majority  of  the  electrons  in  the  bulk  plasma  are  pushed  with  a  larger  implicit  time  step. 
We  will  now  review  the  direct  implicit  particle  simulation  and  describe  the  multi-time  scale  scheme 
in  the  code. 


2  Direct  Implicit  Particle  Simulation 

Direct  implicit  particle  simulation  will  now  be  presented  for  the  case  of  a  one  dimensional  unmag¬ 
netized  electrostatic  plasma.  Following  Langdon,  Cohen,  and  Friedman  [5],  implicit  particle  movers 
advance  the  particle  position  by  the  equation 

i»+i  =  /3Af*a"+‘  (1) 

where  z  is  the  portion  of  the  position  advance  dependent  on  quantities  known  at  time  level  n  and  /9 
depends  on  the  particular  implicit  scheme.  Now  let  the  electric  field  E  be  defined  on  a  spatial  grid. 
Then, 

a"+‘  =  g£;"+»(z’‘+»)/m,  (2) 

where  weighting,  e.g.  linear,  NGP,  etc.,  would  be  used  to  determine  E  at  particle  location  z""''* 
from  E  on  the  grid.  Combining  Eq.  (1)  and  (2) 

z"+*  =  a'r"+‘(z"+*)  -I-  z"+*,  (3) 

where  a'  =  ^At^q/m. 

Note  that  a  logistical  problem  with  the  implicit  method  is  that  depends  on  which  in 
turn  depends  on  the  particle  locations  z"'*'*.  Unfortunately,  all  of  the  particle  locations  depend  on 
One  way  to  get  around  this  problem  is  to  linearize  the  locations  z"'*’’  about  the  locations 
i"+*  as  in  Reference  1.  Then  p"'*'*  on  the  grid  is  written  as 


p»+i 

where  depends  on  and 

ip"+‘  =  -a,[^+‘^z"+‘]. 


(4) 

(5) 


Equation  (5)  is  derived  from  linearizing  the  particle  to  grid  weighting  equations  about  z’''*'^  Strictly 
speaking,  iz"'*'*  =  z""*"*  —  is  an  individual  quantity  for  each  particle  and  is  given  by  Eq.  (3). 
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Instead,  is  approximated  as  a  grid  quantity  while  maintaining  the  form  as  given  by  Eq.  (3). 

That  is, 

(6) 

Equations  (5)  and  (6)  may  be  combined  and  then  substituted  into  the  right  hand  side  of  the  Poisson 
equation  to  get  an  implicit  version 

a.[l  +  'Ao,  (7) 

where  o  =  o'/fo-  The  p  terms  are  determined  by  weighing  the  x  terms  to  the  grid. 

3  Implicit  Particle  Advance,  Field  Solve,  and  Boundary 
Conditions  for  RF  Discharge  Modeling  in  PDPl 

There  are  several  possible  implicit  finite  difference  approximations  to  the  particle  equations  of  mo¬ 
tion.  Here,  Friedman’s  adjustable  damping  scheme  is  chosen  [8].  The  scheme  is 

X-+*  =  i"  -b  (8) 

At [a"+*  +  tfa"/2  -I-  (1  -  «/2)A"-*]/2,  (9) 

A"-’  =  (1  -  «/2)a"**  -I-  ^"-3/2,  (10) 

where  a"’"'*’*  =  is  the  acceleration  of  the  particle  which  is  found  by  interpolating 

the  field  known  on  a  spatial  grid  to  the  particle  location.  The  A  terms  are  the  lag  accelerations 
which  damp  high  frequency  oscillations.  9  =  0  gives  no  damping  and  9=1  gives  the  D1  scheme. 
For  this  scheme,  0  =  1/2. 

Equation  (7),  the  time  implicit  Poisson  equation,  was  written  for  a  single  species.  For  multiple 
species,  Eq.  (7)  is  generalized  to 

5»[1  +  X]dt4>  =  -p/fo.  (11) 

where  the  n  +  1  superscript  is  suppressed  and 

pj  —  E,qEijWj(i),  (12) 

Xj  =  E,(9"AtV2m€o)Ei.yW,(x).  (13) 

As  in  the  previous  section,  x  is  the  portion  of  the  advance  dependent  on  quantities  known  at  time 
level  n.  The  W  terms  indicate  weighing  particles  to  the  grid.  The  finite  difference  version  of  Eq. 
(11)  is 

(1  +  -  [2  +  Xj-t/2  +  +  [1  +  Xj+i/2]4>j+i  =  -Ax^pjfeo,  (14) 
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where  X>±i/2  =  iXj  +  Xi±i)/2i  j  =  1,2,  ...,nc-  2,nc-  1,  and  nc  is  the  number  of  grid  cells.  For 
the  RF  discharge,  left  plate  is  biased  at  the  RF  source  voltage  and  the  right  plate  is  referenced  to  0 
potential.  This  gives  =  Vrf(0  and  <l>j=nc  =  0- 

The  electric  field  at  the  interior  grid  points  can  be  determined  by  the  finite  difference  version  of 
E  =  -dj4),  i.e., 

Ei  =  (<t>j-x-4>i+x)l(2£.x).  (15) 

At  j  =  0  and  j  =  nc,  an  implicit  Gauss  law  must  be  used  to  determine  the  field.  An  implicit  Poisson 
equation  can  be  written  in  vector  notation  as 


V-(l  +  x)E  =  p/fo.  (16) 

The  integral  representation  of  Eq.  (16)  gives  the  implicit  Gauss  law; 

J^(l  +  x)E-dS  =  Q/co,  (17) 

where  Q,  the  approximation  to  the  enclosed  charge,  is  due  to  p  and  the  surface  charge  density. 
Drawing  a  Gaussian  pillbox  for  Eq.  (17)  around  the  left  wall  and  spanning  the  distance  j  =  -1/2 
to  j  =  1/2  gives 

[(1  +  x)-E]j=i/2  =  +  ^Ax/2(o,  (18) 

where  <To  is  the  wall  surface  charge  density.  Inside  the  perfectly  conducting  wall,  £j’<o  =  0.  At  the 
wall.  Gauss  law  gives  Ejso  =  <J'o/<o-  Substituting  this  expression  into  Eq.  (18)  and  solving  for  Ej-o 
gives 

Eq  =  [(1  +  Xi/2)(^o  -  ^i)/Ai]  —  ^Az/2eo,  (19) 

and  similarily,  for  the  right  wall, 


Enc  —  [(1  +  Xnc-l/2)(4inc— I  ^ne)/Al]  +  PncAz/2eo.  (20) 

Finally,  it  is  assumed  that  if  z  penetrates  the  wall,  then  the  particle’s  final  location  will  be  in  the 
wall. 


4  Numerical  Procedure  for  Direct  Implicit 

The  numerical  procedure  for  direct  implicit  simulation  without  multi-time  scale  scheme  is  as  follows. 

•  Take  the  particle  quantities  z",  r",  a"~*,  A"~^  and  the  grid  quantities  E",  p"  as  known. 

•  In  the  prepush  for  each  species,  calculate  A”~^  using  Eq.  (10),  then  calculate  a"  using  linear 
weighting  from  the  grid  to  the  particles,  and  finally,  calculate  z’’"''*  using  the  known  quantities 
in  Eqs.  (8)  and  (9). 
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•  Use  linear  weighting  from  the  particles  to  the  grid  to  determine  and  in  Eqs.  (12) 

and  (13).  Solve  Eq.  (14)  for  ^"+1  with  a  tridiagonal  matrix  inverter.  Calculate  using 
Eqs.  (15),  (19),  and  (20). 

•  In  the  postpush  for  each  species,  determine  and  v""*"*  using  the  terms  in  Eqs.  (8) 
and  (9)  for  all  particles  with  i"'*’*  not  in  the  walls.  Here,  a  gE"‘''^(x"'''*)/m  and  this 
quantity  is  not  saved. 

•  Repack  the  particle  arrays  to  remove  particles  lost  to  the  walls.  Do  Monte  Carlo  collisions. 


5  Multi-Time  Scale  Criteria 

We  now  describe  the  scheme  for  the  electrons  only,  knowing  that  ions  can  also  be  pushed  with  the 
large  implicit  time  step.  The  v,  —  x  phase  space  is  divided  into  two  regions  by  placing  a  rectangle  in 
phase  space.  Inside  the  rectangle,  electrons  are  pushed  with  and  outside  the  rectangle,  electrons 
are  pushed  with  6t  where  At  >  jf.  The  rectangle  boundaries  in  x  are  the  sheath  plasma  boundaries 
averaged  over  an  RF  cycle.  The  rectangle  boundaries  in  v  are  ±£Ax/At  where  0  <  c  <  1.  The 
boundary  in  x  is  necessary  because  the  electric  field  in  the  sheath  has  large  spatial  gradient. 


6  Multi-Time  Scale  Method 

Here,  we  construct  a  very  simple  multiscale  method  for  electrons  as  a  correction  to  our  original 
implicit  scheme.  The  correction  is  an  attempt  to  accurately  model  the  fast  electrons.  Suppose 
At  =  46t,  where  6t  is  chosen  to  accurately  model  the  fast  electrons.  The  fast  electrons  are  advanced 
in  the  time  steps 

n  =  0, 1, 2, 3, 4, 5, 6, 7, 8, ... 

The  slow  electrons  are  advanced  in  the  time  steps 

AT  =  0,4,8,... 

Both  fast  and  slow  particles  are  moved  with  the  direct  implicit  simulation  technique  for  consis¬ 
tency.  If  the  fast  and  slow  electrons  were  moved  independently,  the  implicit  corrections  to  the  free 
stream  densities  are 

Sf,}.  =  -d4{0q/m)^f.6t^Er], 

for  the  fast  (F)  electrons  and 

6p^  =  -dA{0q/m)p^At^E^], 
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for  the  slow  (S)  electrons.  /?  depends  on  the  particular  implicit  scheme  used  to  push  the  electrons. 
The  field  must  be  calculated  at  every  time  step  with  the  Poisson  equation: 

+  ipf- +  ^  +  6p5)/fo- 

The  key  is  to  determine  ps  and  6ps  at  the  n  time  steps  during  which  the  slow  electrons  are  NOT 
pushed,  e.g.  4  to  5,  5  to  6,  etc.  Since  we  view  fast  electron  motion  as  a  correction  to  the  implicit 
scheme,  6ps  retains  the  same  form: 

f>p's  =  -3,[(/39/m)p5Af^r"], 

where  we  keep  the  same  At  as  if  we  were  actually  pushing  the  slow  particles  from  t  —  At  to  t. 

The  approximation  to  ps  for  the  time  steps  in  which  the  slow  electrons  are  not  pushed  is  given 
by  linear  weighting: 

p|  =  (3/4)^  +  (l/4)p|, 

^  =  (2/4)^ +  (2/4)^, 

Ps  =  (1/4)4 +(3/4)^. 

Note  that  4’*  lo  fh*  actual  free  streaming  of  the  slow  electrons.  The  implicit  Poisson 

equation  becomes 

where  x  =  Boundary  conditions  may  be  worked  out  from  an  implicit  Gauss 

law  based  on  the  above  equation. 

7  Numerical  Procedure  for  Multi-Time  Scale 

The  numerical  procedure  for  moving  the  simulation  from  time  step  4  to  8  is  outlined  below.  The 
procedure  is  the  same  for  moving  the  simulation  from  8  to  12,  12  to  16,  etc. 

•  Determine  which  electrons  are  slow  and  which  electrons  bk  fast. 

•  Prepush  slow  electrons  from  4  to  8. 

•  Calculate  ps  for  5,6,7. 

•  Prepush  fast  electrons  from  4  to  5. 


•  Field  solve  at  5. 


•  Postpiish  fast  electrons  from  4  to  5. 

•  Repeat  fast  electron  steps  for  5  to  6,  6  to  7,  7  to  8. 

•  Postpush  slow  electron  from  4  to  8. 
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Fig.  1.  A  parallel  plate  RF  discharge 


Fig.  2.  A  typical  electron  energy  distribution 
function  in  an  RF  discharge 
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Abstract 

Two  methods  are  presented  for  WpcAt  >>  1  particle  simulation.  The  first  method  pertains 
to  cases  where  there  is  no  VteAt/L  <  1  restriction.  For  this  case,  an  alternate  form  of  the 
direct  implicit  particle  simulation  method  is  presented.  The  second  method  pertains  to  cases 
where  there  is  a  VttAi/L  <  1  restriction,  where  L  is  a  sheath  length  that  must  be  resolved  by 
the  particles.  The  method  uses  quasineutrality  in  the  plasma  bulk  with  a  logical  sheath  on  the 
plasma  boundary.  Quasineutrality  removes  the  <  1  restriction  and  the  logical  sheath 

allows  some  resolution  of  the  sheath  physics. 


1  Introduction 

Very  often,  one  wants  to  simulate  a  plasma  phenomenon  with  a  characteristic  frequency  u)  «  Upe- 
However,  if  the  model  equations  include  Upt  physics  as  well  as  w  physics,  such  as  the  equations 
of  electrostatic  particle  simulation,  then  the  Upe  time  scale  must  be  resolved  or  there  will  be  a 
numerical  instability.  Furthermore,  there  are  often  spatial  scales  L  such  as  sheath  lengths  or  gradient 
lengths  that  also  need  to  be  resolved.  For  the  case  of  particle  simulation,  this  gives  the  restriction 
VuAt/L  <  1.  This  restriction  has  nothing  to  do  with  numerical  stability  but  everything  to  do  with 
modeling  the  particular  physical  problem  accurately. 

To  overcome  the  stability  restriction  WpeAf  <<  1,  one  can  either  alter  the  equations  used  to 
model  the  physical  phenomenon  or  use  implicit  methods  to  solve  the  original  equations.  Both 
methods  may  fail  if  there  are  other  restrictions.  For  example,  one  may  introduce  quasineutrality 
into  the  model  to  remove  the  Wp,  time  scale,  but  this  approximation  breaks  down  in  the  sheath  region 
where  the  plasma  is  very  non-neutral.  As  another  example,  one  may  introduce  implicit  methods, 
but  VteAt/L  <  1  must  be  obeyed  for  sheath  resolution  and  if  L  ^  Xpe,  one  still  has  the  UpeAt 
restriction. 
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Here,  two  methods  are  presented  to  allow  large  Wp^At.  The  first  method  is  just  an  alternative 
form  of  direct  implicit  particle  simulation  [?],  but  this  alternative  form  may  be  more  accurate.  This 
method  would  apply  to  the  case  of  an  RF  driven  sheath  simulation.  For  this  problem,  the  plasma 
sheath  region  is  large  due  to  the  Child  Langmuir  law  and  Vte^t/L  <  1  is  not  too  limiting.  The 
second  method  attempts  to  remove  both  the  Upe  time  scale  and  the  L  space  scale  restrictions  by 
using  a  quasineutral  Poisson  equation  [?]  to  model  the  plasma  bulk  and  a  logical  sheath  [?]  to 
model  the  plasma  boundary  layers.  This  method  would  apply  to  the  case  of  a  plasma  between  two 
floating  conductors.  For  this  problem,  the  sheath  is  on  the  order  of  a  few  A/je- 

The  plan  of  this  report  is  as  follows.  In  the  second  section,  an  alternative  form  of  direct  implicit 
particle  simulation  is  presented.  In  the  third  section,  a  numerical  proceedure  for  thb  method,  in 
the  remainder  of  this  report  called  method  a,  is  outlined.  In  the  fourth  section,  the  derivation  of  a 
quasineutral  Poisson  equation  is  given  and  the  highest  frequency  that  the  quasineutral  model  can 
resolve  is  calculated.  In  the  fifth  section,  the  logical  sheath  is  reviewed  and  it  is  shown  how  to  merge 
the  bulk’s  quasineutral  model  with  the  boundary’s  logical  sheath.  In  the  sixth  section,  a  numerical 
proceedure  for  this  method,  in  the  remainder  of  this  report  called  method  b,  is  outlined.  Finally,  in 
the  seventh  section,  some  concluding  remarks  are  made. 


2  Method  A:  An  Alternative  Form  of  Direct  Implicit  Par¬ 
ticle  Simulation 

A  form  of  direct  implicit  particle  simulation  will  be  presented  for  the  case  of  a  one  dimensional 
unmagnetized  electrostatic  plasma.  Following  Langdon,  Cohen,  and  Friedman  [?],  implicit  particle 
movers  advance  the  particle  position  by  the  equation 

i"+i  =  -t-i"'*'’,  (1) 

where  z  is  the  portion  of  the  position  advance  dependent  on  quantitie.s  known  at  time  level  n  and 
depends  on  the  particular  implicit  scheme.  Now  let  the  electric  field  E  be  defined  on  a  spatial  grid. 
Then, 

a"+' =9f7’+‘(x"+*)/m,  (2) 

where  weighting,  e.g.  linear,  NGP,  etc.,  would  be  used  to  determine  E  at  particle  location  z"+* 
from  E  on  the  grid.  Combining  Eq.  (1)  and  (2) 

z"+*  =  o'E"+>(z"+*)  +  z"+',  (3) 


where  o'  —  0At^qfm. 
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Note  that  a  logistical  problem  with  the  implicit  method  is  that  depends  on  p""*"’  which  in 

turn  depends  on  the  particle  locations  x"'*'*.  Unfortunately,  all  of  the  particle  locations  depend  on 
One  way  to  get  around  this  problem  is  to  linearize  the  locations  x""*"*  about  the  locations 
jn+i  ^  jjj  Reference  1.  Then  p""*"*  on  the  grid  is  written  as 

p"+»  =p"+i+ip"+i,  (4) 

where  depends  on  x"'*’*  and 

ip"+i  =  (5) 

Equation  (5)  constitutes  the  linearization.  Strictly  speaking,  ix"+*  =  x"+’  —  x"+*  is  an  individual 
quantity  for  each  particle  and  is  given  by  Eq.  (3).  Instead,  ix”"*"*  is  approximated  as  a  grid  qucintity 
while  maintaining  the  form  as  given  by  Eq.  (3).  That  is, 

6x"+*  =  a'E’'+\  (6) 

Equations  (5)  and  (6)  ,iay  be  combined  and  then  substituted  into  the  right  hand  side  of  the  Poisson 
equation  to  get  an  implicit  version 

a,[i  +  =  -^+V*o,  (7) 

where  a  =  a'/eo. 

Another  approach  will  now  be  described  which  does  not  depend  on  linearization  or  the  approx¬ 
imation  used  to  arrive  at  Eq.  (6).  Suppose  a  one  dimensional  grid  is  defined  with  grid  points  0  to 
N.  Let  /»"+*  be  the  nearest  left  grid  point  for  a  particle  located  at  x"+*,  i.e., 

=  int(x"+V^ar),  (8) 


where  Ax  is  the  space  between  grid  points.  Given  A"'*'*  and  on  the  grid,  Eq.  (3),  with  linear 
weighting,  would  be 


Equation  (9)  can  be  solved  for  VUft.+i^j(x"'’'*)  =  (x""*"*  —  X(,.+i)/Ax,  which  is  the  linear  weighting 
factor  from  a  particle  at  x"+>  to  the  A""''*  +  1  grid  point.  Equation  (9)  becomes 


where  y  =  a/ Ax.  Also, 


be;.*.', + 


VUft.+,{x"+‘)=l-lVfc.+.+,(x"+»)- 


Now  for  any  species,  density  on  the  grid,  is  formed  as 


(12) 
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where  the  sum  is  over  particles  that  contribute  density  to  the  j  grid  point  and  Q  =  q/Ax  for  one 
dimensional  weighting.  Particles  in  the  grid  cells  immediantly  to  the  left  and  right  of  the  grid  point 
j  contribute  the  density.  contributes  from  the  left  grid  cell  and  contributes 

from  the  right  grid  cell.  With  this  in  mind,  Eq.  (12)  can  be  written 

+  (13) 

where  the  first  summation  is  for  contributions  from  the  left  grid  cell  and  the  second  summation  is 
for  contributions  from  the  right  grid  ceil. 


Substituting  Eqs.  (10)  and  (11)  into  Eq.  (13)  gives 


„n+l  _ 


For  the  first  summation,  j  —  1  is  substituted  for  inside  the  summation.  For  the  second  sum¬ 
mation,  j  is  substituted  for  h’'*^  inside  the  summation.  The  subscript  j  and  j  ±  I  field  terms  are 
constant  over  the  summations.  Therefore, 


™  (,  he::,',  + 


(1  -  -  i^.V)]  ) 


+  Q 


where  is  the  number  of  particles  weighing  in  from  the  left  and  is  the  number 

of  particles  weighing  in  from  the  right. 


To  use  Eq.  (15)  to  form  the  density,  one  must  know  before  knowing  i"'*'*.  A  possible 
solution  to  this  problem  is  to  guess  =  tn((x’''*'*/Az).  This  is  a  good  guess  since  particle 
Courant  conditions  require  VtAt/Ax  <  1,  and  so  one  may  expect  |  —  i"'*'*  |<<  Ax.  For 

particles  very  close  to  a  grid  point,  there  may  be  a  jump  from  one  grid  cell  to  the  next  making  the 
above  guess  for  /i"+‘  incorrect.  However,  particles  very  close  to  a  grid  point  are  moved  mostly  by 
a  held  located  at  that  grid  point  and  these  particles  contribute  charge  mostly  to  that  grid  point. 
Therefore,  any  errors  due  to  the  wrong  being  used  only  appear  in  side  terms,  i.e.  a  small  push 
from  and  a  small  contribution  to  j  ^  1  instead  of  j  ±  1.  Linear  weighting,  on  which  this  method  is 
based,  minimizes  these  kinds  of  errors. 


Eq.  (15)  is  substituted  into  the  right  hand  side  of  the  Poisson  equation.  The  resulting  set  of 
equations  would  have  the  form 


d2^"+‘  =  E.p'*+‘(£7+»,r;*+‘), 


(16) 
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E"+‘ = (17) 

From  Eq.  (15),  it  can  be  seen  that  the  particle  quantities  need  only  be  summed  over  once.  Then  Eqs. 
(16)  and  (17)  must  be  solved  simultaneously  over  the  grid.  This  cannot  be  done  by  simple  matrix 
inversion  techniques  because  of  the  nonlinear  form  of  Eq.  (15).  Instead,  some  matrix  iteration 
technique  may  be  used.  There  are  continuing  advances  in  matrix  iteration  techniques  which  may  be 
applied  to  this  case. 

Method  A  has  been  presented  for  a  one  dimensional  system  with  line£ir  weighting  from  the 
particle  to  the  grid  and  visa  versa.  It  is  not  difficult  to  see  how  the  method  may  be  generalized  to 
higher  dimensions,  however,  the  equations  would  be  quite  messy  even  for  two  dimensions.  Using 
higher  order  weighting  would  make  even  the  one  dimesional  equations  cumbersome  also.  Also,  it  is 
difficult  to  see  how  to  incorporate  NGP  weighting  due  to  the  form  of  the  equations  presented.' 

Finally,  one  can  compare  this  method  to  Appendix  C  in  Langdon,  Cohen,  and  Friedman  [?].  In 
that  appendix,  the  authors  present  a  correction  to  the  perturbed  particle  position  which  has  a  form 
similar  to  Eq.  (10).  However,  that  perturbed  particle  position  is  still  substituted  into  the  linearized 
theory  as  outlined  in  Eqs.  (5)-(7).  Here,  an  attempt  is  made  to  use  n  +  1  particle  positions,  written 
in  terms  of  n  +  1  fields  on  a  spatial  grid,  without  appealing  to  a  linearization  form  as  represented  by 
Eq.  (5).  Whether  this  is  just  nitpicking  or  if  this  method  is  worth  the  trouble  remains  to  be  seen. 


3  Method  A:  Numerical  Proceedure 


A  numerical  proceedure  by  which  the  above  method  may  be  implemented  will  now  be  presented. 

For  a  one  dimensional  unmagnetized  electron  and  ion  plasma,  the  following  equations  describe  the 

evolution  of  the  particle  orbits  and  the  electrostatic  field: 

i  =  t>, 

(18) 

V  =  qE{x)/m, 

(19) 

=  -E,p(z)/fo, 

(20) 

E{x)  -  dr<^(x). 

(21) 

The  time  implicit  discretized  form  of  these  equations  may  be  written 

i"+>  =  i"  +  A<(o"+>  +r")/2. 

(22) 

=  r"  +  q£^t[E"+\x^+^)  +  E^(x^)]/2m, 

(23) 

=  -E.p;+‘Ao, 

(24) 
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£:;+»  =  (25) 

where  Eqs.  (22)  and  (23)  are  advanced  for  all  particles  in  both  species  and  Eq.  (15)  is  substituted 
into  the  right  hand  side  of  Eq.  (24).  The  trapezoidal  scheme  for  the  particle  advance  is  chosen  for 
illustrative  purposes  only.  There  are  other  schemes  that  contain  adjustable  damping  [?]. 

The  first  step  of  the  particle  advance  is  to  calculate  the  part  of  the  advance  dependent  on  known 
time  level  n  quantities: 


ti"+‘  =  «"  +  qAtE^{x’')/2Tn, 

(26) 

jn+l  _  J.n  ^  At(ti'*+‘  +  u")/2. 

(27) 

h"+*  ss  inf(f’‘+7Ax). 

(28) 

Next,  the  weighting  factors  in  Eq.  (15)  are  calculated.  Note  that  this  means  two  grid  quantities 
must  be  found.  Particles  weighing  in  from  a  left  grid  cell  to  a  grid  point,  i.e.  the  summation  over 
/,n+i  .f  1  =  is  calculated  as  a  grid  quantity.  Particles  weighing  in  from  a  right  grid  cell  to  a  grid 
point,  i.e.  the  summation  over  =  j,  is  also  calculated  as  a  grid  quantity.  Then  Eq.  (24),  with 
Eq.  (15)  used  on  the  right  hand  side,  and  Eq.  (25)  are  iterated  on  the  grid  to  find  the  potential  and 
the  electric  field. 

Finally,  the  particles  must  complete  their  advance: 

^n+i  jn+i  ^  ,AfE;+‘(i"+‘)/2m,  (29) 

i"+‘  ss  +  Af(»"+*  -  i;"+*)/2.  (30) 

The  scheme  presented  can  be  considered  a  first  step  in  an  interation  over  particle  and  grid  quantities 
if  one  wants  to  make  the  computational  effort.  The  boundary  conditions  can  be  treated  in  the  same 
manner  as  in  explicit  methods. 


4  Method  B:  The  Quasineutral  Model  for  the  Plasma  Bulk 

Now  consider  a  problem  of  a  quasineutral  plasma  bulk.  One  can  expect  important  physics  in  long 
wavelength  low  frequency  modes.  Therefore,  the  plasma  approximation,  n*  =  n;  but  electrostatic 
E  ^  0,  can  be  used  explicitly  in  constructing  a  physical  model.  Hewett  derived  the  quasineutral 
equations  for  a  particle  ion,  fluid  electron,  electromagnetic  Darwin  model  [?].  A  quasineutral  model, 
complete  with  a  quasineutral  Poisson  equation,  will  now  be  derived  for  the  simpler  unmagnetized, 
electrostatic  case  where  both  species  are  treated  as  particles. 
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The  fluid  equations,  in  vector  form,  for  an  unmagnetized,  electrostatic,  non  resistive  plasma  with 
electrons  and  singly  charged  ions  are 


-  edtTig  +  V  .  J«  =  0, 

(31) 

c3,n,  +  V  •  J,  =  0, 

(32) 

dtJg  -  eV  ■  tig  <  vv  >«=  e^ngE/mg, 

(33) 

dtJi  +  cV  •  Ri  <  »»  >,=  e^niE/rrii, 

(34) 

where  J  =  ±enu  and  u  and  <  vv  >  are  first  and  second  moments  of  the  distribution  which,  like  n, 
can  be  found  by  particle  weighting. 

The  basic  idea  behind  the  quasineutral  Poisson  equation  is  to  use  the  plasma  approximation 
rit  =  Rt  directly  and  ignore  the  regular  Poisson  equation  which,  with  this  approximation,  would  give 
=  0.  With  the  plasma  approximation,  the  sum  of  the  continuity  equations,  Eqs.  (31)  and  (32), 

gives 

V.(J.  +  J,)  =  0.  (35) 

With  the  plasma  approximation,  the  sum  of  the  current  equations,  Eqs.  (33)  and  (34),  gives 

fli(Je  +  J.)  +  eV  •  n,(<  vv  >,  -  <vv  >*)  =  e*(l/me  +  l/m,)n,E.  (36) 

With  the  definitions  P  =  V  •  nj(<  vv  >i  -  <  vv  >«)  and  /i  =  mgmi/{Tne  +  m^),  Eq.  (36)  becomes 

a,(Je  +  J0  +  eP  =  c*n.E//i.  (37) 

Operating  on  Eq.  (37)  with  V-  and  cancelling  the  resulting  J  terms  using  Eq.  (35)  gives,  after  some 
manipulation,  the  quasineutral  Poisson  equation: 

=  E  •  Vnj/nj  — /iV  •  P/cRj,  (38) 

where 

E  =  -V^.  (39) 

A  numerical  solution  of  Eqs.  (38)  and  (39)  will  not  require  iteration  on  the  grid  because  the  right 
hand  side  of  Ek].  (38)  is  linear  in  E. 

It  should  be  noted  that  because  the  ions  and  electrons  are  treated  as  particles,  one  could  calculate 
both  Tie  and  rii  on  the  grid.  The  above  theory  only  requires  rii  because  it  is  assumed  that  Re  =  Rj. 
In  a  particle  simulation,  n^  »  Ui  even  in  a  quasineutral  model  because  of  particle  noise.  However, 
in  the  spirit  of  a  low  frequency  long  wavelength  quasineutral  simulation,  n,  as  determined  by  the 
ions  weighed  to  the  grid  will  be  used  in  Eq.  (38)  and  the  definition  of  P.  n,  could  be  calculated 
and  then  compared  to  n^  for  diagnostic  purposes. 
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To  calculate  the  natural  frequency  present  in  the  model  equations,  one  starts  with  Eqs.  (31)  to 
(34),  takes  n*  =  rij,  and  then  perturbs  the  equations  such  that  n  =  rjo  +  and  u  =  6u.  Ignoring 
the  66,  letting  6  oc  exp[i(k  •  x  -  wt)],  and  taking  V  •  n  <  wti  >=  v^Vn,  one  gets  a  linear  dispersion 
relation  of  the  form 

=  k^(T,+Ti/m,  +  mi),  (40) 

which  predicts  ion  acoustic  waves.  These  waves  have  the  highest  frequency  that  needs  to  be  resolved 
by  the  numerical  analogues  of  the  particle  advance  and  quasineutral  Poisson  equations. 


5  Method  B:  The  Logical  Sheath  Model  for  the  Plasma 
Boundary  Layer 

The  quasineutrality  of  the  plasma  ends  near  a  bounding  wall  since  the  more  mobile  electrons  charge 
up  the  wall  while  leaving  the  plasma  near  the  wall  ion  rich.  This  region  is  called  the  sheath  and  is 
usually  a  few  wide  for  an  undriven  wall.  The  sheath  is  analagous  to  boundary  layers  in  fluid 
flow.  In  this  region  quasneutrality  breaks  down.  Furthermore,  even  in  steady  state,  the  sheath 
represents  a  sharp  gradient  and  the  accuracy  condition  Vt^t/L  <  1  forces  small  At  despite  the  large 
At  that  can  be  used  in  the  bulk. 

Fortunately,  the  logical  sheath  boundary  condition  can  be  used  so  that  the  basic  features  of 
the  sheath  can  be  modeled  without  creating  an  accuracy  condition  on  At  [?].  The  logical  sheath 
derivation  is  based  on  the  fact  that  a  potential  drop  near  a  floating  wall  always  adjusts  itself  so  that 
the  electron  and  ion  fluxes  into  the  plate  are  equal  at  steady  state.  The  model  collapses  the  sheath 
potential  drop  into  a  step  function  with  the  sheath  width  approaching  zero.  To  the  plasma  particles 
Just  touching  the  wall,  the  electric  fleld  and  the  wall  charge  still  appear  to  be  zero.  This  allows  the 
boundary  conditions  for  Eqs.  (38)  and  (39)  to  be 
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(41) 

4>o  =  ^ii 

(42) 

where  subscript  0  indicates  the  zeroeth  grid  point  at  the  left  wall.  There  are  similar  conditions  at 
the  right  wail. 

The  electric  fleld  tangential  to  the  wall  is  always  zero  while  the  logical  sheath  allows  the  electric 
field  normal  to  the  wail  to  be  zero  thus  giving  Eq.  (41).  Ek|uation  (42)  proceeds  from  a  Gaussian  pill 
box  drawn  about  the  wall.  The  charge  enclosed  in  that  pill  box  is  zero  because  the  logical  sheath 
suggests  that  there  b  no  surface  charge  density  on  the  wall  and  quasineutrality  suggests  that  there 
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is  no  charge  density  near  wall  or  any  other  place  in  the  plasma.  The  bulk  is  allowed  to  press  up 
against  the  wail. 

In  reality,  there  is  an  electric  field  near  the  wall,  but  only  the  particles  that  can  cross  the  potential 
barrier  and  can  go  into  the  wall  feel  the  field.  This  effect  is  modeled  as  follows.  The  number  of 
electrons  and  the  number  of  ions  Ni  crossing  into  the  wall  are  counted.  If  >  Ni,  the  Ng 
electrons  are  ordered  in  velocity  from  fastest  to  slowest.  All  Ni  ions  and  the  fastest  A,  electrons  are 
absorbed  into  the  wall  to  give  equal  fluxes  while  the  slower  Ng  —  Ni  electrons  are  reflected  as  if  they 
were  in  a  potential  well.  For  the  rare  case  Ng  <  Ni  all  particles  are  absorbed  and  a  positive  wall 
surface  charge  density  is  calculate  which  leads  to  a  modification  in  Eq.  (42): 

<l>o  =  <t>i +cru,Ax/€o-  (43) 

No  modification  in  Eq.  (41)  is  needed  since  the  activity  of  the  electric  field  at  0  is  modeled  by  the 
particles  themselves. 


6  Method  B:  Numerical  Proceedure 

A  numerical  proceedure  for  method  b  will  now  be  presented.  It  will  contain  less  detail  than  the 
proceedure  presented  for  method  a  because  method  b  is  rooted  in  physics  whereas  method  a  is  rooted 
in  numerics.  Furthermore,  the  numerical  proceedure  to  be  presented  is  explicit.  However,  the  logical 
sheath  boundary  condition  has  been  used  with  implicit  codes  that  solve  the  regular  Poisson  equation 
in  the  bulk  [?]  and  there  is  no  reason  why  the  particle  advance  and  quasineutral  Poisson  equations 
cannot  be  solved  by  implicit  numerical  methods. 

Equations  of  motion  for  electrons  and  ions  take  the  usual  form 

X  =  V,  (44) 

V  =  gE(a:)/m,  (45) 

where  E  is  known  on  a  grid  and  the  field  is  interpolated  to  the  particle.  After  the  advance,  the 
logical  sheath  boundary  conditions  are  calculated.  Then,  the  <  vv  >  and  n  terms  are  accumulated 
on  the  grid  and  P  is  calculated.  Eqs.  (38)  and  (39)  can  be  solved  by  matrix  inversions  or  iterative 
methods  so  that  the  field  at  the  next  time  step  is  determined. 


7  Conclusion 

Two  methods  have  been  presented  to  simulate  plasmas  with  the  UpgA  constraint  relaoced.  Method  a 
pertains  to  the  case  where  there  are  no  gradient  constraints  whereas  method  b  pertains  to  the  case 
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of  a  thin  sheath  layer.  Both  methods  would  require  ViAt/Ax  <  0.5  for  accuracy.  It  now  rem2Lins  for 
both  methods  to  be  tried  for  some  simple  bounded  plasma  test  cases. 
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A  Comparison  of  PIC  Simulation  and  Experimental  Results  in  a 

Capacitative  RF  Discharge  [1] 

P.  Mirraishidi,  B.P.  Wood,  V.  Vahedi,  G.  DiPeso 
September  24,  1991 

Simulation  results  from  PDPl,  a  ld3v  bounded  particle-in-cell  code[2],  are  compared  to  recently 
published  experimental  results[3]  over  a  pressure  range  of  lO-lOOmTorr  and  lOO-lOOOV  applied  RF 
voltage  in  a  symmetric,  parallel  plate,  Argon  discharge.  We  show  that  where  similar  results  are  ob¬ 
tained,  the  simulation  allows  insight  into  plasma  parameters  which  are  not  experimentally  accessible, 
such  as  details  of  the  electron  power  sources  and  losses. 

We  have  been  exploring  the  numerical  effects  of  various  algorithms  of  moving  particles.  ,  In 
particular,  how  employing  certain  algorithms  lead  to  either  numerical  cooling  or  heating.  Implicit 
codes,  for  example,  are  not  able  to  model  the  collision  of  the  fast  tail  of  electrons  in  a  RF  discharge 
with  the  plasma  sheath  accurately.  Over  long  periods  of  time,  the  lack  of  stochastic  heating  in  the 
simulated  plasma  will  result  in  a  completely  inaccurate  system.  The  large  time  steps  associated  with 
implicit  codes  curtails  their  use  in  the  simulation  of  RF  discharges. 

Below  are  two  examples  of  phase-space  snapshots  of  PDPl,  showing  both  an  explicit  result(physically 
correct)  and  an  implicit  result(which  is  numerically  cooled)  for  RF  discharges;  note  that  the  latter 
does  not  have  any  fast  electrons. 
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Progress  on  PDP2,  A  Two  Dimensional  Bounded  Particle 

Simulation  Code 


Gregory  DiPeso  and  Vahid  Vahedi 
September  30,  1991 


Abstract 

PDP2,  a  two  dimensional  bounded  electrostatic  particle  simulation  code  for  plasmas,  is  a 
fairly  straightforward  extension  of  the  one  dimensional  bounded  particle  simulation  code  PDPl. 
In  this  report,  we  will  describe  in  detail  the  extension  of  the  field  and  circuit  solve  to  a  second 
dimension  in  space  with  periodic  boundary  conditions.  We  will  also  discuss  particle  loading, 
pushing,  and  weighting  in  two  dimensions.  Like  the  current  workstation  version  of  PDPl,  PDP2 
is  run  in  the  X windows  environment  and  so  we  will  describe  the  available  diagnostics.  Finally, 
we  will  conclude  with  the  additional  physics  and  diagnostics  that  should  be  added  before  PDP2 
can  acheive  its  full  potential. 


1  Introduction 

Two  dimensional  particle  simulation  of  plasmas  is  not  new.  Previous  examples  of  two  dimensional 
particle  models  in  the  electrostatic  limit  include  the  work  of  Decyk  and  Dawson  [1]  and  the  work  of 
Thielhaber  and  Birdsall  [2].  Decyk  and  Dawson  developed  a  Poisson  solver  for  generalized  boundary 
conditions  based  on  two  dimensional  Fourier  transforms.  They  successfully  simulated  surface  waves 
at  a  vacuum  plasma  interface.  Thielhaber  and  Birdsall  developed  a  Poisson  solver  for  a  plasma 
bounded  by  perfectly  conducting  walls  in  x  and  periodic  in  y  [3]  in  their  code  ES2  to  study  Kelvin 
Helmholtz  modes. 

The  model  used  in  PDP2  is  basically  the  same  as  that  used  in  ES2.  The  differences  between 
PDP2  and  ES2  are  inherited  from  PDPl.  The  physics  differences  are  that  PDP2  will  use  the  charged 
particle  neutral  collision  model  of  PDPl  [4]  and  uses  the  external  circuit  solver  of  PDPl  [5]  where 
as  ES2  uses  collisions  only  as  a  source  of  electrons  and  ions  to  counter  particle  losses  to  the  walls 
and  uses  the  external  circuit  solver  of  Lawson  [6].  The  PDPl  collisional  model  is  an  improvement  in 
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that  it  treats  collisions  in  a  more  sophisticated  manner.  The  PDPl  circuit  solver  is  an  improvement 
in  that  it  allows  a  greater  range  of  input  parameters.  The  diagnostic  differences  are  that  PDP2  runs 
in  an  Xwindow  environment  like  PDPl  [8].  This  allows  the  viewing  of  physics  as  it  happens  and 
interactive  measurements  of  any  diagnostic  via  crossh?=>-  and  rescaling  mechanisms.  ES2  runs  on 
the  CRAY  in  which  graphics  files  can  only  be  viewed  at  the  end  of  execution. 

The  development  of  PDP2  is  incomplete.  The  actual  addition  of  the  collisional  package  and  the 
corresponding  changes  to  input  parameter  reads  are  required.  The  options  of  particle  injection  from 
the  walls  and  distributed  sources  could  be  added.  PDPl’s  dumping  and  restarting  capability  should 
be  added.  Furthermore,  all  of  the  PDPl  diagnostics  must  be  included  in  PDP2.  Nonetheless,  PDP2 
has  already  been  used  to  sucessfully  simulate  sheath  waves  in  an  unmagnetized  plasma  [9].  Perhaps 
a  more  obvious  test  would  be  to  temporarily  use  the  ES2  collision  package  and  try  to  recover' the 
results  of  Thielhaber  and  Biidsall. 


2  Field  and  Circuit  Solve 

Consider  a  plasma  bounded  in  x  by  perfectly  conducting  walls  and  assume  the  plasma  is  periodic  in 
y  as  in  Fig.  1.  The  conducting  walls  are  attached  to  an  RLC  external  circuit  with  a  voltage  source. 
The  wall  to  wall  distance  is  Lf  Let  grid  numbers  run  from  i  =  0,  M  in  x  and  j  =  0,N  in  y.  The 
walls  are  at  grid  point  i  =  0  and  i  —  M.  Zero  net  Charge  over  the  entire  system  gives  the  relation 
N-l  N-i  M 

'^{AyL.clj  +  AyL,<T*„^)+  ^  ^AxAyL,p\j  =  0,  (1) 

J=0  >=0  tiiO 

where  the  superscript  denotes  the  time  level.  The  surface  charge  density  a  varies  in  y  because  of  its 
relation  to  the  internal  charge  density  p.  Since  the  z  dimension  does  not  figure  into  the  model,  L, 
may  be  chosen  to  be  the  unit  length,  e.g.  Im  in  MKS.  Next,  Kirchhoff’s  Laws  are  employed  at  the 

left  t  =  0  wall  and  the  right  i  =  M  wall  which  gives 

N-\  N-l 

qU  +  Qco  -  Q‘o^‘  =  E  -  E  (2) 

j=0  >=0 

and 

N-l  N-l 

QpM  +  =  E  -  E  (3) 

}=0  >=0 

where  Qp  is  the  net  charge  extracted  or  deposited  on  the  wall  due  to  absorption  of  plasma  charged 
particles  or  emission  of  charge  particles  to  the  plasma  and  Qc  is  the  charge  on  the  external  capacitor. 
Note  that  Qeu  is  available  from  Qco  by  adding  Eqs.  (2)  and  (3)  and  substituting  Eq.  (1). 

The  external  circuit  equation  is 


I'QcO  +  RQcO  +  Qeo/C  =  V(f)  +  +  (f>0i 
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where  <p  indicates  an  equipotential  at  the  wall  (no  j  subscript)  and  F  is  a  voltage  source.  The 
solution  to  this  equation  can  be  written  as  [5] 

Q[o  =  {Vit)-K*-\-<l>*M-4>o)/oo,  (5) 

where  ao  and  K*  are  numerical  factors  defined  in  [5]. 

To  develop  boundary  conditions  in  x  for  the  Poisson  equation,  we  draw  a  Gaussian  pillbox  about 
some  grid  point  j.  This  is  shown  on  Fig.  2  for  the  left  i  =  0  wall.  Noting  that  the  electric  field 
directed  along  the  wall  (along  y)  and  inside  the  wall  is  zero  and  that  at  i  =  1/2  is  given  by 
~  finite  differenced  version  of  Gauss’  Law  is 

(K  -  =  Aa;(/9o,i/2  +  <Tj,j/Ax)/eo,  (6) 

where  the  superscript  indicating  the  time  level  is  dropped  for  internal  4)  and  p  .  For  the  right  hand 
wall,  a  similar  equation  may  be  written  as 

(<^M  -  <^M-i,i)/Ax  =  Ax(pmj72  +  ffMj/Ax)/fo.  (7) 

The  boundary  conditions  for  y  are  periodic.  Finally,  the  Poisson  equation  is  written  numerically  as 

-  24>i,j  +  +  (^i,i-i  -  +  ^i,>+i)/Ay*  =  -pi,}/(o  (8) 

for  the  internal  <f)-  This  completes  the  fundamental  set  of  equations  for  the  held  and  circuit  solve 
routine.  For  more  discussion  see  Chap.  16  in  Birdsall  and  Langdon  [3]  and  Lawson  [7]. 

We  solve  the  set  of  above  coupled  equations  by  Fourier  transform  methods  as  discussed  in  Chap. 
14  of  Birdsall  and  Langdon  [3].  The  y  periodic  boundary  conditions  are  taken  into  account  natu¬ 
rally  and  wave  spectrum  (in  ky)  diagnostics  are  a  by-product  of  the  transformations.  The  Fourier 
transform  of  Eq.  (8)  in  the  y  direction  is 

+  Dk<i>i,k  +  <t>i+l,k  =  -Ax^p,,t/fo,  (9) 

where  the  subscript  k  is  the  same  as  ky  and  =  — 2[l-f2(Axsin(7rm/A^))^]  is  a  term  that  accounts 
for  the  finite  difference  terms  in  Eq.  (8)  [3],  <t>i,k  is  the  Fourier  transform  of  4>ij,  it  >  0,  and 
«  =  1,  M  —  1.  Since  (p  is  constant  in  space  along  the  conducting  walls,  e.g.,  =  4>o  for  all  j,  then 

4>o,k  =  0  for  fc  >  0  and  <t>o,k=o  =  0o-  Similarily,  =  0  for  ik  >  0  and  4>o,k=o  =  is 

generally  a  complex  quantity  where  as  <t>ij  is  real.  The  k  =  0  mode  Fourier  transformed  phi  has  a 
zero  imaginary  part  and  the  real  part  is  just  the  average  over  the  periodic  length  Ly.  The  Fourier 
transform  of  Eq.  (6)  for  <:  =  0  is 


(<^0  ~  <i>\,k=o)  =  Ax^(po.k=o/2  +  (Tq  ii-o/Ax)/fo, 


(10) 
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where 

N-\ 

i=o 

We  do  not  need  the  equations  for  the  higher  modes  of  o. 

Substituting  Eqs.  (5)  and  (11)  into  Eq.  (2)  gives 

<i>o'¥ckQ^yL,N<TQi^_Q  —  0^,  (12) 

where  /3‘  =  V(t)  —  A'‘  +  ato(Qpo  ~  equation  relates  the  average  left 

wall  surface  charge  density  <To,t=o  1°  potential  )/>q.  Rearranging  Eq.  (10)  gives 

^0  ~  ^i,t=o  -  ^a;<To,t=oAo  =  Ai*po,t=o/(2eo)-  (13) 

Finally,  writing  Eq.  (9)  for  ib  =  0  gives 

^•-i,t=o  +  Dk=o<t>i,k=o  +  ^t+i,t=o  =  —  Ax*p,,t=o/fo»  (14) 

where  i  =  l,Af  —  1.  Recall  that  the  ib  =  0  mode  Fourier  transforms  are  real  quantities.  For  i  =  1, 
^o,k=o  =  ^0  as  usual  and  for  i  =  M  —  1,  <i>M,k=Q  =  4>*]^  =  0  as  a  reference  potential.  There  are  M  +  1 
unknowns  which  are  <To,fc=o'  ^O'  ^i.ksOi  *  =  1,M  -  1.  Equations  (12)-(14)  constitute  M  +  1 
equations  for  all  of  the  unknowns.  The  resulting  matrix  for  this  system  is  tridiagonal  which  is  easily 
inverted.  The  above  set  of  equations  couples  the  zero  mode  portion  of  the  Poisson  equation  to  the 
circuit  equation.  For  the  ib  >  0  modes,  no  coupling  to  the  circuit  is  present  since  <f>o^k  =  4>M,k  =  0 
for  k  >  0.  For  this  case,  Eq.  (9)  gives  M  —  1  equations  for  the  M  -  1  unknowns  ib  =  1,  Af  —  1. 
Good  algorithms  for  fast  Fourier  transforms  of  real  quantities  and  tridiagonal  matrix  inversion  can 
be  found  in  Numerical  Recipes  in  C  by  Press  et  al  [10].  These  algorithms  are  used  in  PDP2. 

Finally,  the  can  be  inverse  Fourier  transformed  into  Eqs.  (6)  and  (7)  can  be  then  used 
to  get  ffoj  and 

^oj  -  ~  -  Ax(po,>/2),  (15) 

=  fo('^M  -  -  Az(pm,>/2).  (16) 

Once  (Toj  and  ^M,j  are  calculated,  E,  on  the  wall  is  given  by, 

Ex,o,j  =  o’o.yAoi  =  <^o,>Ao-  (17) 

=  0  on  the  wall,  smd  for  the  internal  points,  the  electric  field  is  obtained  from  the 
finite  difference  version  of  E  =  — V0, 

(Ex)ij  =  —  ^,+i)y/(2Ax),  (18) 


—  (^>-1  ~  ^>+i)«/(2Aj/). 


(19) 
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At  y  =  0  =  Ly,  periodicity  is  used  to  complete  the  finite  difference.  For  a  short  circuit,  4>q  =  V{t) 
and  =  0  as  a  reference  potential.  Then,  Eq.  (14)  can  be  solved  independently  of  Eqs.  (12)  and 
(13).  For  an  open  circuit,  the  Qeo  terms  in  Eq.  (2)  vanish,  leaving 

^o,k=o  —  Qpo/(^y^*^)t  (20) 

so  that  Eqs.  (13)  and  (14)  can  be  solved  with  Eq.  (17)  instead  of  Eq.  (12). 


3  Particle  Manipulation 

A  particle  in  a  simulation  is  often  called  a  superparticle  since  it  represents  many  actual  plasma 
particles.  For  instance,  if  a  plasma  chamber  volume  is  Im^  with  a  density  of  10*®m~^  fo  ach  of 
two  species,  then  there  are  2  x  10**  plasma  particles  in  the  chamber.  For  even  the  fastest  computers 
with  the  largest  storage  space,  it  is  impossible  to  represent  the  actual  number  of  particles.  If 
we  use  200, 000  simulation  particles  to  represent  the  system,  then  each  of  these  superparticles  is 
worth  10**  plasma  particles.  Despite  the  seemingly  poor  representation,  many  theoretical  and 
experimental  results  have  been  observed  via  simulation.  In  this  section,  the  word  particle  will  refer 
to  the  simulation  superparticles. 

Manipulation  of  the  particles  consists  of  loading,  solving  of  the  orbit  equations,  and  calculating 
density  on  the  grid  from  the  particle  positions.  Particle  injection  and  particle  collisions,  not  yet 
implemented  in  PDP2,  are  discussed  elsewhere  [3,  4]. 

Loading  of  particles  consists  of  choosing  the  initial  conditions  for  the  problem  of  interest  and 
then  setting  the  initial  locations  and  velocities  of  all  of  the  particles  so  as  to  represent  the  initial 
state.  In  PDP2,  we  assume  an  initial  drifting  Maxwellian  distribution  in  velocity  and  a  uniform 
distribution  in  space.  The  Maxwellian  takes  the  usual  isotropic  form 

/(»)  =  C  exp(-»*/2v|-),  (21) 

where  v,  is  the  thermal  spread  in  velocity  and  t?*  =  (r,  -  v,o)*  +  (v*  -  v,o)*  +  (Vr  -  v,o)^.  The  zero 
subscripted  values  are  the  user  specified  drift  velocities.  Once  Eq.  (18)  is  inverted  to  yield  v  [3],  the 
velocities  in  each  direction  are  chosen  via  regular  random  numbers  and  the  drifts  are  then  added. 
The  uniform  loading  in  configuration  space  is  done  by  using  bit  reversed  random  numbers  in  y  and 
random  numbers  in  x.  At  this  time,  there  is  no  option  to  add  a  perturbation  to  the  particle  orbits. 
Any  waves  that  are  to  be  excited  must  come  from  the  thermal  noise  in  velocity. 

Once  the  particles  are  loaded,  their  equations  of  motion  must  be  solved  to  advance  the  velocities 
and  positions  in  time  by  one  At.  Then  charge  density  on  the  grid  is  accumulated  (using  linear 
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weighting)  so  that  the  Poisson  equation  may  be  solved  and  the  electric  field  can  be  calculated  for 
the  next  time  advance.  The  equations  of  motion  are  simply 

X  =  V,  (22) 

V  =  {g/m)[E(x)  +  v  x  Bo],  (23) 

where  Bo  is  a  user  specified  magnetic  field.  The  user  may  specify  the  magnitude  and  direction 
of  the  magnetic  field.  Equations  (19)  and  (20)  are  advanced  by  the  well  known  Boris  mover  [3]. 
Particles  that  go  beyond  y  =  0  or  y  =  £y,  are  simply  reintroduced  into  the  other  end  in  y  to 
satisfy  the  periodic  boundary  conditions.  Particles  that  hit  the  wall  on  the  left  side  contribute  to 
Qpo-  Since  Qpu  is  not  used,  the  particles  that  hit  the  wall  on  the  right  side  do  not  contribute  to 
the  physics  except  indirectly  through  the  circuit  equations.  Any  particle  that  hits  a  wall  is  deleted 
from  the  particle  list.  To  reduce  the  number  of  multiplications  and  divisions,  the  following  internal 
normalizations  are  used: 

q*-qNc/{L,^xAy), 

X  x/Ai, 
y  4-  y/Ay, 

V,  4—  VfAt/Ax, 

Vp  4—  «yAt/Ay, 

V,  *-  v„ 

where  Nc  is  the  ratio  of  plasma  particles  to  superparticles. 

Finally,  one  needs  to  calculate  the  charge  density  on  the  grid.  To  do  this,  simple  linear  weighting 
is  used.  It  should  be  mentioned  that  linear  weighting  is  used  to  determine  E  at  the  particle  for  the 
Boris  mover.  From  Fig.  3,  the  fraction  of  the  particle  charge  contributed  to  grid  point  /  (or  the 
fraction  of  electric  field  contributed  from  grid  point  I  as  needed  by  the  Boris  mover)  is  given  by  the 
ratio  of  the  area  of  square  /  to  the  total  cell  area  AxAy  where  1=1,2, 3, 4. 

We  have  also  added  the  option  of  a  drift  kinetic  species.  In  this  approximation,  the  gyromotion 
of  the  superparticles  is  neglected.  The  equations  of  motion  take  the  form 

X  =  Ei/B  +  w||B/B,  (24) 

«||  =  (25) 

where  Ej.  =  E  x  B/B,  =  E  •  B/B,  and  B  is  the  magnitude  of  B.  Note  that  the  electron 
superparticles  follow  motions  parallel  to  B  and  the  E  x  B  drift.  Since  there  are  no  gradients  or 
curves  in  the  magnetic  field,  VB  and  curvature  drifts  are  absent  from  the  equations.  E  varies  in 
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space  which  gives  a  correction  to  the  E  x  B  drift  on  the  order  of  For  small  Vce,  this  correction 
ran  be  ignored.  E  also  varies  in  time  which  sets  up  a  polarization  drift  which  is  not  accounted  for 
in  the  Eqs.  (24)  and  (25). 

The  numerical  method  used  to  advance  Eqs.  (24)  and  (25)  cannot  be  a  simple  leapfrog  because 
the  advance  of  x  depends  on  the  field  quantities  directly.  Here,  we  propose  the  following  multistep 
centered  leapfrog  scheme 


‘  =  x"-*  -t-  +  2El(x’*)]/B, 

(26) 

+  9A<E||(x")/m, 

(27) 

x"-‘  =  x"+‘  Af[B»y+‘/*]/B, 

(28) 

where  Ex  and  E\\  are  defined  on  the  grid.  Note  that  fields  calculated  at  x"  are  used  to  move  both 
t)||  from  n  —  1/2  to  n  +  1/2  and  x  from  n  —  1  to  n  +  1.  In  the  next  pass  through  this  algorithm,  x" 
is  moved  to  x"'*'*,  is  moved  to  and  the  fields  used  in  these  moves  are  calculated  at 

xn+i  algorithm  has  been  programmed  but  has  yet  to  be  tested  for  accuracy  and  stability. 


4  Concluding  Remarks:  A  Wish  List 

The  basic  changes  to  PDP2  tire  in  order  of  importance; 

•  Add  the  charge  particle  neutral  collision  package  and  tdl  of  the  input  deck  and  reading  para- 
phanalia  that  is  needed. 

•  Add  all  of  the  diagnostics  already  present  in  PDPl.  This  includes  number  versus  time  for  any 
species  which  is  a  very  basic  and  important  diagnostic. 

•  Incorporate  the  new  version  of  Xgraphics  including  three  dimensional  visualization. 

•  Add  the  file  dumping  and  restarting  options  already  present  in  PDPl. 

•  Include  particle  injection  from  either  wall. 

Other  changes  to  PDP2  are  conceivable.  For  example,  one  may  want  to  consider  the  case  of  walls 
in  the  y  as  well  as  the  x  direction.  Then  one  could  use  either  the  field  solve  of  Decyk  and  Dawson  or 
opt  for  the  many  elliptic  solvers  currently  avjulable.  If  PDP2  is  to  be  run  on  a  parallel  computing 
machine,  obviously  one  would  want  to  choose  a  Poisson  solver  that  takes  the  most  advantage  of 
parallel  architecture.  Depending  on  the  application,  one  may  want  to  model  walls  that  are  not 
perfect  conductors  or  walls  which  change  material  properties  along  their  length.  A  wall  with  a  time 
and  space  varying  potential  could  be  added.  Cylindrical  and  spherical  shapes  are  possible. 
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PDP2  could  be  made  into  a  direct  implicit  or  even  a  multiscale  code.  The  collisonal  model 
could  be  made  even  more  sophisticated  by  taking  into  to  account  charged  particle  collisions.  The 
possibilities  are  of  course  endless,  but  with  its  present  sophistication,  PDP2  could  be  a  good  base 
from  which  other  particle  simulation  codes  can  be  developed. 
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Theory  and  Simulation  of  Plasma  Sheath  Waves 
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Abstract 

Sheath  waves  have  been  investigated  analyticaily  and  with  particle  simulation  for  an  nn- 
magnetized  two  dimensional  plasma  slab  with  periodic  boundary  conditions  in  y  and  conduct¬ 
ing  walls  in  z.  Analytically  treating  the  sheath  as  a  vacuum  layer,  the  sheath  wave  bears  a 
resemblance  to  plasma  vacuum  surface  waves.  The  simulations  are  in  good  agreement  with  the 
theory  for  both  bulk  Bohm  Gross  waves  and  edge  sheath  waves. 


1  Introduction 

It  is  well  known  that  there  is  a  great  variety  of  waves  in  a  plasma  that  is  well  neutralized  (iti  ~  rit) 
and  does  not  have  sharp  gradients  in  field  or  density  quantities.  Waves  also  exist  at  the  plasma 
edge  or  sheath  where  there  is  large  charge  unbalance  (n,-  jt  Uf)  and  where  the  gradient  scale  lengths 
are  on  the  order  of  the  electron  Debye  length  in  the  unmagnetized  case  or  on  the  order  of  the  ion 
gyroradius  in  the  magnetized  case.  These  waves  have  received  less  attention  in  the  literature.  This 
paper  b  a  report  on  electrostatic  waves  propagating  along  the  unmagnetized  plasma  edge  or  sheath. 
Both  analytic  theory  and  computer  simulation  are  used  to  study  the  sheath  waves.  The  computer 
simulation  may  be  viewed  as  an  experiment  if  the  simulation  model  b  constructed  from  first  principal 
physics  with  a  minimum  of  approximations  or  assumptions. 

Before  we  start  on  the  two  dimensional  theory  and  simulations,  let  us  review  the  results  of  one 
dimensional  simulations  [1].  The  one  dimensional  simulations  are  bounded  by  perfectly  conducting 
wails  which  are  connected  by  an  external  RLC  circuit  with  optional  voltage  and  current  sources. 
The  simplest  boundary  conditions  for  which  sheath  formation  is  observed  b  the  short  circuit  where 
the  reference  potential  or  voltage  on  both  walls  b  set  to  zero.  The  device  is  initially  filled  with 
warm  electrons  at  a  density  neo-  The  electrons  have  a  full  Maxwellian  velocity  dbtribution  at  a 
temperature  T«.  The  ions  are  treated  as  an  immobile  background  with  a  density ’n,o  =  ngo  so  that 
the  system  is  initially  neutral.  The  device  length  is  about  50Ax>e. 
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We  connected  the  2d3v  electrostatic  magnetized  bounded  particle  simulation 
code  ES2B,  which  has  been  developed  at  Tohoku  University  in  Japan,  to  Xgrafix. 
ES2B  is  an  x-y  two  dimensional  code.  A  uniform  magnetic  field  is  pointing  in 
the  positive  x-direction.  The  electrons  and  ions  are  continuously  emitted  from  the 
boundary  at  jc  =jci,  where  xi  is  the  system  length  in  x-direction.  The  potential  at 
X  =x\,y  =0,  and  y  =  y  i  are  zero,  where  y  i  is  the  system  length  in  a  y-direction. 
We  can  arbitrarily  specify  the  potential  at  the  boundary  at  x  =0.  Poisson’s  equa¬ 
tion  is  solved  by  the  method  of  superposition.  We  modified  the  code  so  that  we 
can  apply  a  time  varying  potential  at  the  boundary  at  x  =0.  The  Berkeley  group 
will  now  use  this  code  for  fully  bounded  plasmas  and  as  a  key  to  an  RZ  model. 


PROGRESS  REPORT:  RZ  Project  ^ 

SUBJECT:  Two  Dimensional,  Cylindrical  (in  r  and  z  directions)  Multi  grid  Poisson  Solver 
BY:  Conway  Gee 


Motivation 

The  ultimate  goal  will  be  to  write  a  computer  code  that  will  simulate  plasma  in  two  dimensional 
cylindrical  coordinates.  Many  physical  systems  are  in  cylindrical  coordinates. 

An  interesting  mathematical  problem  has  arisen  from  the  project.  And  the  problem  is  coming  up 
with  a  ‘fast’  method  to  solve  the  two  dimensional,  cylindric^  Poisson’s  equation  in  the  r  and  z 
directions. 


i  M(r.z))  +  d^0(r,z)/8z^  =  -p(r,z),  where  p  is  the  known  charge, 
r  dr  dr  and  0  the  potential  to  be  solved, 

and  zero  value  boundary  conditions. 

The  system  looks  like  this. 


At  first  we  looked  at  a  Fast  Sine  (Fourier)  Transform  to  find  the  solution  along  the  z  direction 
and  perhaps  a  ‘  fast’  Bessel  series  in  the  r  direction.  The  FST  requires  0(N  logN)  arithmetic 
operations,  which  is  nearly  optimal  in  terms  of  efficiency  (i.e.  it  is  very  fast).  The  algorithm  for 
the  FST  is  already  well  known.  However,  the  suggested  ‘fast’  Bessel  series  will  require  some 
study  for  it  is  not  so  well  known. 

The  disadvantage  of  using  the  FST  is  that  it  is  rather  specialized.  The  FST  can  be  applied 
primarily  to  systems  which  arise  from  separable  self-adjoint  boundary  value  problems. 
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The  Multgrid  Melhcxl 


Introduction 

The  mulngnd  method  is  a  culmination  of  results  from  trial  and  error  ever  since  iterative,  or 
relaxation,  methods  were  first  used  to  achieve  a  solution.  Ideas  about  making  the  iterative  methods 
converge  faster  were  known  for  many  years,  but  only  until  recently  has  these  ideas  been 
formalized  as  the  multigrid  method. 

Review  of  the  Classical  Methods 

If  the  onginal  operator  matrix  A  in  the  system  Ax  =  b  is  too  difficult  to  work  with  (via  Gaussian 
elimination,  etc.),  then  another,  simpler  matrix  can  be  used  in  place  of  A.  Calling  the  new  matrix 
M,  we  have 


M\  =(M-A)x  +  b  (1) 

In  the  classical  iterative  methods,  (1)  is  solved  iteratively  by  successive  substitutions. 
Beginning  with  an  initial  guess  Xq,  it  may  be  the  vector  of  zeros,  the  current  guess  Xj.  leads  to  the 

next  approximation. 


A/Xj.^j  =(M -A  )xj.+ b  (2) 

The  basic  choices  of  M  are;  \.  M  =  diagonal  part  of  A  (Jacobi’s  method) 

2.  M  =  triangular  part  of  A  (Gauss-Seidel  method) 

Note  that  the  choices  of  M  are  ‘close’  to  A.  Why  they  should  be  ‘close’  to  A  has  to  do  with 
convergence. 

The  convergence,  or  divergence,  to  the  solution  x  depends  on  M  and  A.  Subtracting  (2)  from 
( 1),  we  have  the  error  equation  for  e^.  =  x  -  Xj.. 

Mcj.^1  =  (M  -A)cy. 
or  e|.^j  =  M  '(M  -/4)ej. 


=  Be^.,  where  B  =  M  '\M  -  A) 

=  I  -  mU 

The  current  error  is  related  to  the  original  error,  ej.  =  B*^eQ  . 

Thus,  the  goal  is  to  make  ej.  converge  to  zero,  and  that  will  happen  only  when  the  powers  of  B 
converge  to  zero. 


ej. ->  0  and  Xj. ->  x  iff  /?*'->  0 

Recall  that  B  ^Cq  splits  into  terms  of  the  eigenvectors  and  eigenvalues  of  B.  Then 

B'^Cq  -  Cj(lj)*^Wj  +  .  .  .  +  Cj^(lj^)’"Wjj ,  where  1  denotes  the  eigenvalues 

and  w  the  eigenvectors. 


Therefore,  the  powers  converge  to  zero  if  and  only  if  every  cigcin  aluc  of  B  satisfies  the 
condition 


I1J<  1. 

The  rate  on  convergence  depends  on  the  largest  I  Ij  I,  which  is  called  the  spectral  radius  of  B. 
The  matrix  M  should  satisfy  two  criteria: 

1.  The  equation  M  =  (A/  -  /1)X|.  +  b  should  be  easy  to  solve.  Hence,  M  should  be 
convenient  to  use. 

2.  M  should  be  ‘close’  to  A,  so  that  the  eigenvalues  of  B  =  M'\M  -  A)  =  1  -  M  ‘  *  A  are  small 
(less  than  one). 

Jacobi’s  method  takes  Mas  the  diagonal  part  of  A.  If  all  ajj  are  nonzero,  then  the  diagonal 
matrix,  D,  is  easily  invertible.  The  matrix  representaion  is 

=(D  -  A)X|j  +  b. 


ar-d  by  elements 


aii(Xi)k+l  =  (-ai2X2  -  2,3X3 ....  -  ai„x„),.  +  b. 

Wk+1  =  -  ^02X2  -  •  •  •  -  ^nn-l  Vl)k  +  • 

The  Gauss-Seidel  method  is  different  in  that  it  starts  using  each  component  of  the  new  x,.^,  as 
soon  as  it  is  computed.  Then  x^^^,  replaces  x,.  an  element  at  a  time,  and  the  old  vector  \  is 
replaced  as  fast  as  Xj^^,  is  created.  The  first  step  finds  the  first  component  as  in  Jacobi's  method, 
and  the  next  step  operates  immediately  with  that  first  new  component. 

^2(^2)k+l  =  '^I^^l^k+l  *  •  •  •  '  ^n'n^k  ^2 

The  last  equation  uses  the  new  values  exclusively. 

=  (At^l  ■  ^112^2  -  ■  •  •  -  ann-lXn-l)k+l  + 

M  here  is  the  lower  triangular  part  of  A,  when  all  the  terms  x,^^,  are  moved  to  the  left  side.  On  the 

right  side,  M  -  A  is  a  strictly  upper  triangular  matrix. 

Improvements  of  these  methods,  such  as  successive  overrelaxation  and  rcd-black  ordering,  are 
jilso  well  known. 

The  problem  with  the  classical  iterative  methods  mentioned  abo\e,  even  with  improvements,  is 

that  convergence  is  controlled  by  the  largest  eigenvalue  of  B  =  M  "'(M  -  A),  which  usually  is 
associated  with  errors  of  lowest  frequencies.  High  frequency  errors  are  quickly  eliminated,  but  the 
smooth  component  of  error  holds  back  convergence. 

As  an  example  of  the  relationship  between  error  and  eigenvalues,  let's  examine  the  one 
dimensional  scrond  order  differential  equation 

u”(x)  =  -f(x),  0<x<  1,  with  the  boundary  conditions  u(())=u(  1  )=() 


The  domain  of  the  problem  is  partitioned  into  N  subinierv'als  by  the  grid  points 


Xj  =  jh,  where  h  =  1/N  is  the  constant  width  of  the  subintcrv  als.  The  discretized  form  of  the 
equation  is 

(Uj.j  +  2Uj  +  Uj^j)/h2=  -f(Xj) 

Recall  that  the  Jacobi  method  uses  the  diagonal  part  of  the  operator  matrix  such  that 
B  =  I  -DU 

and  written  out  for  our  example,  the  matrix  D  looks  like  the  follow  ing. 


The  eigenvalues  of  B  and  A  are  related  by 
l(fl )  =  1  -  KA )  /  2. 

The  eigenvalues  of  A  are  found  to  be 

1(A  )  =  4sin^(lcni:  /  2N),  where  1  ^  k  £  N-1, 
and  the  eigenvectors,  represented  as  Wj^j,  the  jih  component  of  the  kih  eigenvector, 

w,.j  =  sin(jkK  /  N),  1  ^  k  ^  N-1,  0  ^  j  ^  N. 

Thus,  the  eigenvalues  of  B  are 

ljj(B  )  =  1  -  2sin^(lcrt  /  2N),  1  <  k  <  N-1. 

Modes  in  the  lower  half  of  the  spectrum,  where  the  wave  number  has  the  range  1  s  k  <  N/2,  are 
said  to  be  the  low  frequency  or  smooth  modes.  The  modes  in  the  upper  half  of  the  spectrum,  with 
N/2  ^  k  ^  N-1,  are  called  high  frequency  modes. 

For  the  lowest  frequency,  k  =  1 , 

Ij  =  1  -  2sin^(n  /  2N)  =  1  -  2sin^(7th  /  2)  ~  1  -  rr^h^  /  2, 

and  recalling  the  relationship, 

ej.  =  =  Cjfljl'^w  I  +  .  . .  +  Cjj(l^)*^w here,  k  is  not  the  wave  number 

the  equations  suggests  that  the  eigenvalue  associated  with  the  smcxilhest  mcxle  w  ill  be  very  clo.se  to 
one,  and  therefore  convergence  to  the  solution  w  ill  be  slowed  dow  n.  Notice  that  the  smiler  the 
grid  spacing  h,  the  closer 


I j  is  to  1.  Hence,  any  attempt  to  improve  the  accuracy  by  decreasing  the  si/c  of  the  grid  spacing 
will  only  worsen  the  convergence  of  the  smooth  components  of  error. 

Multi  grid 

The  multigrid  method  changes  the  scale  of  the  problem  so  that  the  smooth  component  of  error 
can  be  reduced  more  rapidly.  Thus,  the  idea  is  to  change  the  grid.  Smooth  errors  on  a  grid  of 
width  h  can  be  attacked  on  a  coarser  grid,  say,  of  width  2h,  where  the  error  is  not  so  smooth.  The 
following  is  an  outline  of  an  elementary  multi  grid  method.  Here,  Q*'  denotes  the  grid  associated 
with  spacing  h,  and  is  the  vector  associated  with  Q^. 

Coarse  Grid  Correction  Scheme,  <-  CG(x*’,b*’) 

Relax  }4  times  on  /i**x**  =  on  with  any  initial  guess. 

(j4  is  any  arbitrary  small  number,  i.e.  3, 4,  or  5) 

Compute  the  residual  r*'  =  b  -  A^x^. 

Transfer  this  vector  to  a  coarse  grid,  r^  =  Cr**. 

Solve  the  coarse  system  =  r^^  on 

Transfer  the  correction  to  the  fine  grid,  e^  =  Fe^*'. 

Make  the  correction  to  x,  =  x^  +  e^. 

Relax  on  A^x^  =  b**  on  until  desired  solution  is  reached. 

The  coarsening  matrix  C  combines  the  nine  values  on  the  fine  grid  to  give  a  single  value,  at  their 
center,  on  the  coarse  grid.  The  weights  are  shown  in  the  figure  and  are  mutiplied  by  1/16.  The 
refining  matrix  F  reverses  what  C  does;  each  coarse  value  splits  into  nine  values  on  the  fine  grid. 


The  next  trick  is  evident.  Why  not  nest  the  CG  scheme?  The  coarse  system  _  ^2h 

also  be  coarsened  further,  ‘a  correction  on  the  correction,’  so  to  speak.  This  recursive  nesting 
makes  the  method  possible  for  asymptotic  optimization.  The  nested  scheme  is  sketched  in  the 
following  diagram,  and  by  its  shape,  it  is  called  the  V-cycle. 

h 

21, 

4/1 

8/i 

lf)li 


A  more  sophisticated  scheme,  called  the  full  multigrid  scheme,  nests  if  V-cycles  themselves. 

The  scheme  begins  by  initializing  b^,  b^, . . and  sets  x^,  x^**, ...  to  zero.  The  following  is  the 
schedule  of  the  full  multigrid  scheme. 
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The  Poisson  Solver 

In  our  case,  the  discretized  del  square  operator  of  the  Poisson’s  equation  has  the  form 
40(1  jVAr^-  (4/Ai^  +  2/Az2)0(O,j)  +  (0(O,j+l)+0(O.j-l))/Az2  =  -poj,forr  =  0, 

(0(ij+l)+0(ij-l))/Az2  +  (£±^0(i+lj)  +  fr-Ar/2)0(i-l.i)  -  {2/Ai^  +  2/Az2)0(ij) 

rAi^  rAi^ 


=  -p(i,j).  forr?sO, 

which  determines  the  operator  matrix.  The  grid  that  describes  the  charges  p  is  initially  coarsened 
so  that  the  full  multigrid  scheme  can  be  used.  The  operator  matrix  remains  the  same  (i.e.  its 
elements  are  not  weighted). 

The  advantage  of  the  iterative  methods,  even  without  multigrid,  is  that  the  boundary  conditions 
can  vary.  Say  the  walls  of  the  cylinder  are  grounded,  or  zero,  then  the  0  matrix  will  have  its 
appropriate  elements  fixed  to  zero.  Even  these  elements  of  zeros  can  be  fixed  to  any  arbitrary  value 
other  than  zero.  Different  geometries  can  ^llso  be  easily  done  with  any  number  of  boundary 
conditions.  All  that  is  needed  to  be  done  is  to  fix  the  correct  0  elements  with  the  proper  v^ues. 

Included  with  this  report  is  the  early  version  of  my  program  using  the  full  multigrid  scheme, 
and  initial  results  appear  promising.  In  one  test,  the  time  it  takes  full  multigrid  to  solve  a  system  up 
to  five  significant  digits  is  roughly  the  same  as  the  time  it  lakes  for  the  FST  to  solve  a  one¬ 
dimensional  problem  with  the  same  number  of  grid  points  (256x256=65536).  In  another  test,  the 

method  solves  for  the  4x4  cylindrical  system  p  =  4(z^  -  Z]Z)/rj^  +  2(r^/r|^  -  1),  where  the  solution  is 

0=(t^/r|^-l)(z^-Z|Z).  The  red-black  Gauss-Seidel  method  takes  roughly  90  fine  grid  iterations  to 

achieve  what  full  multigrid  does  in  only  8  fine  grid  iterations,  along  with  the  much  less  expensive 
coarse  grid  iterations. 


A  New  Approach  to  Traveling  Wave  Tube  Simulation  and  Design[l] 

E.H.  Chao,  C.K.  BirdsaU 
September  26,  1991 


IfiC  (Interactive  Beam  Circuit)  is  a  one-dimensional  traveling  wave  tube  simulation  that  runs  on 
PC'S  and  UNIX  workstations  under  X.  The  code  uses  PIC  techniques  to  model  the  motion  of 
beam  electrons  and  simple  finite-difference  methods  to  model  the  fields  of  the  coupled  slow-wave 
transmission  line.[2]  We  are  currently  running  simulations  at  differing  values  of  the  Pierce  variables 
QC,  C,  and  b  and  measuring  the  power  gain  at  each  of  these  conditions.  These  measurements  are 
then  compared  to  values  calculated  by  BirdsaU  and  Brewer.  [3],[4]  The  initial  compeirisons  have 
been  encouraging.  In  addition,  changes  are  being  mewle  to  the  oode  to  make  it  more  valuable  to 
tube  designers.  These  changes  include  adding  new  diagnostics  (f(E)  at  collector),  modifying  old 
diagnostics,  and  optimizing  the  code. 

Shown  below  are  some  diagnostics  for  C  =  0.096,  QC  =  0.0799,  and  b  =  0.521. 
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Abstract 


Heavy  particles  may  play  a  role  in  determining  the  average  potentials 
experienced  by  ions  in  RF  discharges,  hence  ion  acceleration  into  targets. 
Particulates  or  dust  particles  also  may  play  a  role  in  many  other  plasmas. 
Hence,  it  is  desirable  to  fmd  where  these  heavy  particles  reside  (with 
respect  to  Ae  edge  of  the  plasma,  or  sheath)  and  Iheir  effect  on  Ae  time- 
average  potential  which  accelerates  ions  through  the  sheath. 

Using  our  many-particle  PIC-MCC  ld3v  simulation  code  PDPl^^,  we 
have  been  able  to  show  that  the  particulates  tend  to  become  charged  nega¬ 
tively,  using  cross  sections  for  electron  and  ion  attachment  worked  out  by 
Trulsen^'®^  inspired  by  work  of  R.  Carlile  at  Univ.  Of  Arizona^^l  We  have 
placed  one  heavy  particle  at  various  locations  in  the  sheath  and  found  the 
location  where  the  time  average  field  is  a  minimum;  this  is  then  the  resi¬ 
dence  of  the  particle,  which  turns  out  to  be  very  near  the  time  averaged 
sheath  edge.  We  are  now  putting  in  a  large  number  of  particulates  and 
allowing  them  to  reach  their  respective  equilibrium  positions.  Ion  trans¬ 
port  is  observed  only  after  the  addition  of  a  viscous  term  in  the  force  equa¬ 
tion.  Without  the  addition  of  the  ion  drag,  the  dust  particles  exhibit 
behaviors  similar  to  those  of  the  negative  ions.  However,  with  the  viscous 
term  in  the  force  equation,  the  particulates  exhibit  behaviors  consistent 
with  those  observed  in  the  laboratory. 


1  Current  address:  Department  of  Physics,  UCLA,  Los  Angeles,  CA  90024. 

2  Current  Address:  Univ.  of  Tromso,  P.O.  Box  953,  N9001,  Tromso,  Norway. 
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Introduction 


•  A  dust  particle  will  tend  to  charge  negatively  and  therefore  acquire  an  electric 
potential  negative  charge.  Intuitively,  the  excess  charge  is  caused  by  the 
higher  thermal  velocity  of  the  electrons. 


Due  to  its  excess  negative  charge,  particulates  wiU  create  a  local  potential  bar¬ 
rier  inside  the  plasma.  This  will  cause  many  interesting  effects: 

•  Scattering  of  waves  with  wavelength  X<Xp 

•  Effects  the  overall  charge  balance  of  the  system. 

•  Creates  a  potential  barrier  for  positive  ions  (See  Fig.  1). 

•  Competes  with  other  collisional  processes  in  a  laboratory  plasma. 


Our  interest  in  the  dusty  plasma  is  related  to  the  presence  of  sub-micron  dust 
contamination  in  plasma  aided  processing.  The  dust  particles  are  produced  by 
chemical  and  mechanical  means.  In  some  cases,  nucleation  and  growth  from 
plasma  negative  ions  and  etch  products  is  the  mechanism  for  particle  forma¬ 
tion.  In  o^er  cases,  stress-inducing  processes  may  fracture  thin  films  on 
chamber  surfaces  thereby  injecting  particles  into  the  plasma.  Laboratory 
studies  show  that  these  particles  are  suspended  near  the  sheath  edge  and  drop 
onto  the  wafer  when  the  rf  power  is  turned  off,  thereby  contaminating  critical 
product  surfaces.  In  the  semiconductor  industry,  it  has  been  claimed  that  up 
to  50%  of  rejection  rale  is  directly  caused  by  the  bombardment  of  dust  par¬ 
ticles  on  the  wafers^^^  Our  main  objective  is  to  develop  a  self-consistent 
model  to  simulate  dust  particles  in  RF-driven  discharges.  In  addition,  we  will 
show  that  our  model  agrees  with  basic  probe  theories  when  applied  to  undri¬ 
ven  discharges. 
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PIC-MCC  Scheme 


From  optical  measurements  made  by  Singh  et.  we  have  gained  some 
fundamental  understanding  of  the  particulates.  In  the  situations  which  we  are 
interested,  it  has  been  observed  that  the  particulates  (dust  particles)  are: 

•  Spherical  in  shape. 

•  Lx)cated  at  the  sheath  edge. 


So,  given  that  the  dusty  particles  are  spherical  in  shape,  we  can  defme  the 
capturing  cross  section  as  the  cross-section  for  absorption  of  electrons  and 
ions  by  Ae  dust  particle  (particulate),  i.e.: 


This  model  has  been  implemented  using  a  Paiticle-in-Cell,  Monte-Carlo  col¬ 
lision  (PIC-MCC)  scheme,  i  a  Monte-Carlo  scheme,  the  probability  for 
each  particle  to  undergo  collision  in  a  time-step  A/  is  given  by: 


-e 


In  our  code,  we  generate  a  random  number  (r)  between  0  and  1 ,  then,  the  ran¬ 
dom  number  is  compared  with  the  collisional  probability(P(v)).  If  the  ran¬ 
dom  number  falls  within  the  collisional  probability,  then  a  collision  is  to  take 
place  within  that  time-step.  In  the  case  of  electron  and  ion  absorptions,  two 
events  take  place  when  a  collision  occurs: 

•  The  charge  of  the  dust  particle  is  incremented  by  an  amount  of  q, 
(charge  of  &e  absorbed  species.) 

•  The  particle  is  annihilated.  (See  Fig.  3) 
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The  advantages  of  this  algorithm  are: 


•  Can  je  easOy  converted  into  2-D  and  3-D  particle  simulations. 

•  Can  be  added  into  an  existing  code  without  changing  the  structure  of  the 
code. 


However,  this  algorithm  also  has  its  flaws: 


•  Each  dust  particle  is  treated  as  a  separate  collisional  process,  hence,  for 
multiple  dust  particles,  the  simulation  can  be  time  consuming. 

•  The  code  requires  s  large  number  of  electrons  and  ions  in  order  to  ^lave 
good  statistics. 
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Simulation  Parameters 


The  collisional  model  has  been  implemented  using  an  one-dimensional,  bounded 
simulation  code,  PDPl.  At  t  =  0,  a  maxwellianized  plasma  uniformly  fills  the 
space  with  temperature  Tq  and  density  n^.  A  sheet  of  dust  particle  is  immersed  in 
the  plasma  at  distance  L^.  The  plasma  is  allowed  to  reach  equilibrium.  The  com¬ 
mon  parameters  are: 

•  Argon  Plasma 

•  Temperature  (at  t=0)  =  leV  to  5eV 

•  Initial  Plasma  Density  =  10’^  m'^ 

•  Neutral  Density  =  0 

•  Lengths  10  cm 

•  Number  of  Simulation  Particles  (t=0)  =  8(XX)  per  species 

In  the  undriven  simulations,  we  set  the  neutral  pressure  to  zero  in  order  to  suppress 
all  competing  collision  processes.  Therefore,  only  attachment  of  electrons  and  ions 
is  allowed  during  these  runs. 

In  the  RF-driven  runs,  neutral  gas  is  turned  on  so  electron-neutral,  ion-neutral,  and 
dust  attachment  are  present  during  each  simulation.  Additionally,  a  500-Volt  RF 
source  is  applied  at  13.56MHz.  TTiis  creates  a  sheath  of  1.70  to  2.60  centimeters, 
depending  on  the  phase  of  the  RF  cycle.  The  additional  parameters  are: 

•  Neutral  Pressure  =  10  mTorr 

•  Frequency  =  13.56  MHz  =  8.2  •  lO’sec’’) 

•  RF  voltage  amplitude  =  500  Volts 
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Resuts: 


PART  I:  Undriven  Discharges 


In  order  to  verify  our  collision  scheme,  we  have  placed  a  sheet  of  dust  par¬ 
ticles  in  an  undriven,  maxwellianized  plasma,  then  allow  the  sheet  to  reach 
equilibrium.  The  floating  potential  of  the  dust  is  then  compared  with  probe 
theory  to  verify  our  algorithm.  The  floating  potential  of  a  probe  can  be  calcu¬ 
lated  as  followed.  Let: 

•  V,  =  Velocity  of  each  species. 

•  T,  =  Temperature  of  each  species. 

In  an  unmagnetized  plasma,  the  plasma  species  have  no  directional  prefer¬ 
ence,  hence,  each  of  the  plasma  species  follows  an  isotropic  Boltzman  distri¬ 
bution,  i.e.: 
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So,  using  the  distribution  above,  we  can  calculate  the  current  for  the  plasma 
species: 


nv)=mv) 

$ 


I,iV)  =  nqA 


dv,e 


X 

■5*r 


Here  we  assumed  the  probe  is  facing  the  x-direction,  hence,  only  the 
x-component  of  the  velocity  contributes  to  the  current.  However,  for  elec¬ 
trons,  &e  above  equations  must  be  modified,  because  the  slow  electrons  are 
repelled  by  the  probes.  So,  if  we  defme: 
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The  quantity  v^„  is  the  threshold  velocity  for  the  electrons  and  the  slow  elec¬ 
trons  do  not  contribute  to  the  current.  Therefore,  the  modified  equation  for 
electron  current  becomes: 


J,{V)  =  nqA  r 

Jymin 


2^ 

mv 

_ y 


The  above  description  of  electron  and  ion  behavior  would  yield  the  familiar 
curve  of  the  langmuir  probe.  For  our  purpose,  we  can  think  of  the  dust  par¬ 
ticles  as  floating  probes.  And  the  dust  particle  acquires  a  negative  potential 
such  that  the  equilibrium  electron  flux  equals  to  the  ion  flux. 


IEEE-91 


So,  we  can  derive  the  floating  potential  by  setting  /,  =/,.  This  gives  the  sim¬ 
ple  and  familiar  expression  for  the  floating  potential,  i.e.: 


We  have  compared  the  above  results  with  those  obtained  by  simulation,  and 
we’ve  found  our  model  shows  reasonable  agreement  with  theory. 


Temperature  [eV] 

Vj  (Theory)  [V] 

Vj  (Simulation)  [V] 

1 

-4.32 

-4.18 

2 

-8.64 

-8.23 

3 

-12.96 

-12.21 

4 

-17.28 

-16.59 

5 

-21.60 

-20.85 

The  discrepancies  in  the  results  indicates  that  there  are  excess  electrons  in  our 
model.  However,  if  we  taken  the  spherical  geometry  of  the  probe  into 
account  and  define  the  enhanced  probe  area  through  the  following  equations: 


( 


^i»r.(V,V)=7W' 


,  2eV 

1+ - =rr 

m\v  f 


\ 


( 


A/yy)-%a‘ 


leW  " 


V 


2ey  " 


Hence,  using  the  improved  model,  we  can  recalculate  the  floating  potential, 
once  again,  by  setting  /,  =  /„  and  we  have: 


Temperature  [eV] 

V,  (Theory)  [V] 

Vf  (Simulation)  [V] 

1 

-4.18 

-4.18 

2 

-8.25 

-8.23 
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3 

4 

5 


-12.35 

-16.95 

-20.90 


-12.21 

-16.59 

-20-85 


So,  taken  the  geometry  of  the  probe  into  account,  we  have  shown  that  our 
model  agrees  well  with  theory.  (See  below.) 


Floating  Potential  Measurements 


Planar  Probe  Spherical  Probes  Simulation 
— B —  — —  —-0 
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Part  n:  RF-Driven  Plasmas 


So,  given  that  our  algorithm  is  consistent  with  basic  probe  theories,  we  now  turn 
our  attention  to  RF-driven  discharges.  In  these  mns,  we  are  particularly  interested 
in  two  aspects  of  the  dust  dynamics: 

•  The  position  near  the  sheath-plasma  boundary  which  the  local  time- 
averaged  field  is  a  minimum. 

•  The  macroscopic  dynamics  of  dust  particles  inside  the  bulk. 

To  understand  the  equilibrium  position  of  the  particulates,  we  use  techniques  simi¬ 
lar  to  those  used  in  the  undriven  plasma.  A  sheet  of  immobile  dust  particle  is 
placed  near  the  sheath  of  the  plasma.  Diagnostics  are  added  to  measure  the 
peak-to-peak  average  field  at  the  position  of  the  dust.  The  parameters  are  the  same 
to  those  in  the  undriven  plasma,  with  the  addition  of  a  500  Volt  RF  source  oscillat¬ 
ing  at  13.56MHz.  This  RF  source  produces  a  sheath  ranging  from  1.7  to  2.6  centi¬ 
meters,  depending  on  the  phase  of  the  RF  cycle. 

In  these  runs,  a  large  local  field  is  created  by  the  excess  negative  charges  of  the 
dust  particles.  This  field  creates  a  potential  barrier  hence  effects  the  dynamics  of 
ions  near  the  wall  (i.e.  wafer).  The  large  field  also  changes  the  field  configuration 
around  the  sheath.  So,  by  placing  the  dust  particles  at  various  positions  in  the  RF 
sheath,  we  can  find  the  special  position  where  the  average  field  is  at  a  minimum, 
and  that  is  the  equilibrium  position  of  the  particulates.  Our  measurement  yields  an 
equilibrium  position  of  2.55cm,  which  is  at  the  sheath-plasma  boundary,  consistent 
with  observations  made  in  the  laboratory. 


To  study  the  macroscopic  dynamics  of  the  dust  particles,  we  made  runs  with  a 
large  number  of  dust  particles  in  the  bulk.  The  dust  particles  are  give  finite  mass, 
given  by  the  equation: 
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where  p  is  the  density  of  the  dust  particles  (we  use  water  density  here),  and  is 
the  radius  of  the  dust  particles.  Using  optical  measurements  made  in  the  labora¬ 
tory,  we  choose  the  dust  radius  to  be  0.2  microns.  This  gives  a  mass  of  8  •  10”^' 

grams,  which  is  10^  to  10^  greater  than  the  ion  mass.  The  equilibrium  charge  of  the 
dust  particles  are  also  in  the  order  of  l(X)s  of  the  electrons  and  the  argon  ions, 
hence,  the  dust  charge  to  mass  ratio  of  the  dust  particle  is  in  the  same  order  as  that 
of  the  argon  ions.  In  these  runs,  the  dust  particles  exhibit  behaviors  similar  to  those 
of  the  negative  ions  and  they  are  evenly  distributed  inside  the  bulk  even  after 
hundreds  of  RF  cycles. 

Although  we  have  found  the  equilibrium  position  for  the  particulates  near  the 
sheath,  the  electric  field  around  the  equilibrium  position  is  quite  large  hence  it  is 
very  difficult  for  a  dust  particle  of  fmite  mass  to  maintain  its  equilibrium  position. 
Therefore,  it  is  necessary  to  introduce  additional  forces  that  would  damp  die  the 
large  time-averaged  electric  field  in  the  plasma  sheath.  The  work  done  by  Som- 
merer  et.  al.  suggests  that  ion  drag  force  (including  electrostatic  interactions  and 
interactions  wiA  the  neutral  gas)  plays  a  significant  role  in  the  transport  of  dust 
particles.  In  the  simulation  works  done  by  Sommerer  et  al.,  the  shielded  charge  of 
the  dust  particle  Z’^  is  modified  according  to  the  equation. 


Where  f  =  |  E{z,t)IE^  |,  where  E„  is  the  macroscopic  electric  field  at  a  distance  of 
one  Debye  length  X  from  the  dust  particle  of  charge  Z^.  The  shielded  charge  Z'p 
varies  according  to  the  local  electric  field  and  plasma  conditions.  In  weak  field 
regions  such  as  the  bulk  plasma,  Debye  shielding  is  quite  effective,  and  Z’^  is 

assumed  to  be  1.  Debye  shielding  is  less  effective  for  particles  in  the  strong  RF 
fields  in  the  sheath  because  the  response  times  for  electrons  and  ions  in  the  Debye 
sphere  are  much  longer  than  the  RF  period.  Particles  in  the  sheaths  may  therefore 
respond  to  the  field  as  though  they  are  less  than  fully  unshielded.  The  net  effect  of 
the  modified  Z’^  is  a  strong  confming  force  at  f  ie  sheath  of  the  plasma,  causing 
the  dust  particles  to  stay  near  the  boundary  of  the  sheath. 

In  our  simulations,  we  modify  the  expression  for  the  acceleration  experienced  by 
the  dust  particles  by  adding  a  viscous  term  inside  the  sheath,  i.e.: 


V  -qlm  •£(jc,r)-'n  - V 
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Where  the  viscous  term  "n  includes  effects  such  as  the  interaction  of  dust  particles 
with  the  neutral  gas,  as  well  as  the  electrostatic  interactions  described  in  Sommer- 
ei^^l  Although  the  interaction  of  dust  with  the  neutral  gas  is  unknown,  it  is  intu¬ 
itive  that  the  collision  probability  increases  when  the  particulate  velocity  is  large. 
Inside  the  bulk,  the  velocity  of  a  dust  particle  is  very  small  so  the  particulate  is  not 
likely  to  interact  with  the  neutral  gas.  Therefore,  thermal  motion  dominates  the 
dust  dynamics  inside  the  bulk.  )\^en  the  dust  particle  passes  through  the  sheath,  it 
is  accelerated  and  picks  up  momentum  very  quickly.  We  expect  the  viscous  term 
to  be  significant  when  the  particulate  velocity  is  large,  because  the  dust  particle  is 
more  likely  to  interact  with  the  neutral  gas;  therefore,  the  drag  term  is  included 
only  inside  the  sheath. 

So,  we  add  the  following  lines  to  the  mover  (isp  is  the  species  index  and  i  is  the 
index  for  the  particles): 


/*  PLASMA_SP  is  the  number  of  plasma  species  ♦/ 
/*  LFT_EDGE  and  RHT_EDGE  arc  measured  using 
uncontaminated  RF  discharges  *! 


if((isp<PLASMA  SP)  II 

(x[isp][il>LFT_EDGE  &&  x[ispl[il<RHT_EDGE)) 
v[isp)[i]+=atemp; 

}else 

( 

v[isp][i]+=atemp-cta*v[isp][i]; 


Using  Ti  of  0.01  to  0.03,  we  observe  migration  of  dust  particles  toward  the  edge 
after  approximately  100  RF  cycles. 
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Conclusion 


In  this  paper,  we  presented  simulation  results  for  both  driven  and  undriven  plas¬ 
mas.  In  the  situations  where  there  is  only  one  sheet  present,  it  is  observed  that  the 
particulates  acquire  a  negative  potential  consistent  with  basic  probe  theories. 
Furthermore,  in  the  RF-driven  cases,  we  observe  that  the  average  field  is  at  a  mini¬ 
mum  near  the  edge  of  the  plasma  sheath,  which  agrees  with  experimental  observa¬ 
tions.  Furthermore,  using  water  density,  the  equilibrium  charge-to-mass  ratio  of 
the  dust  particles  is  similar  to  those  of  the  ions;  therefore,  the  dust  particles  would 
behave  like  negative  ions  if  Lorentz  force  is  dominant.  In  the  runs  where  a  large 
number  of  particulates  are  used,  we  believe  that  ion  drag  force,  including  interac¬ 
tion  with  the  neutral  gas  as  weU  as  the  electrostatic  interaction  with  the  plasma  spe¬ 
cies,  is  the  dominant  force  in  dust  transport. 
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Collisions 


Rowchait  of  the  MCC-PIC  scheme  with  the  addition  of  dust  particles 


The  presence  of  the  dust  particles  effects  the  overall  charge  neutrality,  and  creates  a  potential 
barrier  (lower  plot).  The  above  plots  are  made  by  placing  an  immobile  sheet  of  dust  particle  near 
the  sheath  of  the  plasma.  In  this  case,  the  particulate  is  placed  at  =  1.9cm . 
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DnuN^MiMtlr  »c|j) 


The  presence  of  the  dust  particles  effects  the  overall  charge  neutrality,  and  creates  a  potential 
barrier  for  ions  Gower  plot).  The  above  plots  are  made  by  placing  an  immobile  sheet  of  dust 
particle  near  the  sheath  of  die  plasma.  In  this  case,  the  particulate  is  placed  atx^  -  l.Scm. 
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E  Field  vs.  Particulate  Position 


Gas  Pressure  =  10  mTorr 
Frequency  s  13.56  MHZ 


Viscous  force  plays  an  important  role  in  the  transport  of  dust  particles.  The  upper  plot  is  made 
without  drag  force,  and  the  lower  plot  is  the  same  simulation,  with  the  addition  of  a  viscous  term 
in  the  force  equation. 
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